From 98a36a888d2361fcbf4ceb36557a861f263ae248 Mon Sep 17 00:00:00 2001 From: Cedric Muenger Date: Thu, 30 Apr 2026 14:25:37 +0200 Subject: [PATCH 1/6] generated GWP shapefunctions --- src/bases/local/gwpdivlocal.jl | 11 +- src/bases/local/gwplocal.jl | 205 ++++++++++++++++----------------- src/utils/lagpolys.jl | 43 +++++-- 3 files changed, 136 insertions(+), 123 deletions(-) diff --git a/src/bases/local/gwpdivlocal.jl b/src/bases/local/gwpdivlocal.jl index f7d87a82..1c8d024b 100644 --- a/src/bases/local/gwpdivlocal.jl +++ b/src/bases/local/gwpdivlocal.jl @@ -1,6 +1,4 @@ -struct GWPDivRefSpace{T,Degree} <: RefSpace{T} - storage::Vector{@NamedTuple{value::SVector{2,T}, curl::T}} -end +struct GWPDivRefSpace{T,Degree} <: RefSpace{T} end function numfunctions(x::GWPDivRefSpace{<:Any,D}, dom::CompScienceMeshes.ReferenceSimplex{2}) where {D} @@ -11,13 +9,8 @@ function dimtype(x::GWPDivRefSpace{<:Any,D}, Val{(D+1)*(D+3)} end -function GWPDivRefSpace{T,Degree}() where {T,Degree} - NT = @NamedTuple{value::T, curl::SVector{2,T}} - GWPDivRefSpace{T,Degree}(Vector{NT}()) -end - function (ϕ::GWPDivRefSpace{T,Degree})(p) where {T,Degree} - ψ = GWPCurlRefSpace{T,Degree}(ϕ.storage) + ψ = GWPCurlRefSpace{T,Degree}() # (ϕ.storage) dom = domain(chart(p)) u = parametric(p) vals = ψ(dom, u) diff --git a/src/bases/local/gwplocal.jl b/src/bases/local/gwplocal.jl index 9d9cdd1b..296cf378 100644 --- a/src/bases/local/gwplocal.jl +++ b/src/bases/local/gwplocal.jl @@ -1,6 +1,4 @@ -struct GWPCurlRefSpace{T,Degree} <: RefSpace{T} - storage::Vector{@NamedTuple{value::SVector{2,T}, curl::T}} -end +struct GWPCurlRefSpace{T,Degree} <: RefSpace{T} end function numfunctions(x::GWPCurlRefSpace{<:Any,D}, dom::CompScienceMeshes.ReferenceSimplex{2}) where {D} @@ -8,12 +6,7 @@ function numfunctions(x::GWPCurlRefSpace{<:Any,D}, end function dimtype(x::GWPCurlRefSpace{<:Any,D}, dom::CompScienceMeshes.ReferenceSimplex{2}) where {D} - Val{(D+1)*(D+3)} -end - -function GWPCurlRefSpace{T,Degree}() where {T,Degree} - NT = @NamedTuple{value::T, curl::SVector{2,T}} - GWPCurlRefSpace{T,Degree}(Vector{NT}()) + Val((D+1)*(D+3)) end function (ϕ::GWPCurlRefSpace{T,Degree})(p) where {T,Degree} @@ -29,122 +22,122 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di ϕ(dom, u, dimtype(ϕ,dom)) end -function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Dim}, - u, ::Type{Val{NF}}) where {T,Deg,Dim,NF} - - d = Deg - u, v = u - w = 1-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) +@generated function (::GWPCurlRefSpace{T,Degree})(dom::CompScienceMeshes.ReferenceSimplex{Dim}, + bary, ::Val{NF}) where {T,Degree,Dim,NF} + + u = :(bary[1]) + v = :(bary[2]) + w = :(one(T) - $u - $v) + + #nodes for sylvester polynomials + s = :() + for i in 0:Degree+2 + s = :($s..., $i/($Degree+2)) + end - P = SVector{2,T} - NT = @NamedTuple{value::P, curl::T} - # nts = Vector{NT}(undef, NF) - # nts = MVector{NF,NT}(undef) - nts = resize!(ϕ.storage, NF) - idx = 1 + nd1 = :(SVector{2}( -$v, $u-one(T) )) + nd2 = :(SVector{2}( -$v+one(T), $u )) + nd3 = :(SVector{2}( -$v, $u )) + + nts = :() i = 0 - for j in 1:d+1 - k = (d+2)-i-j - Rᵢ = _sylpoly(s, i+1, u) - Rⱼ = _sylpoly_shift(s, j+1, v) - Rₖ = _sylpoly_shift(s, k+1, w) - - 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ₖ - - nts[idx] = (value=Rᵢ*Rⱼ*Rₖ*nd1, curl=curl) - idx += 1 + for j in 1:Degree+1 + k = (Degree+2)-i-j + Rᵢ = gen_sylpoly(s, i+1, u, T) + Rⱼ = gen_sylpoly_shift(s, j+1, v, T) + Rₖ = gen_sylpoly_shift(s, k+1, w, T) + + dRᵢ = gen_sylpoly_diff(s, i+1, u, T) + dRⱼ = gen_sylpoly_shift_diff(s, j+1, v, T) + dRₖ = gen_sylpoly_shift_diff(s, k+1, w, T) + 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ₖ ) + + nts = :($nts..., (value=$Rᵢ*$Rⱼ*$Rₖ*$nd1, curl=$curl)) end - for i in 1:d+1 + for i in 1:Degree+1 j = 0 - k = (d+2)-i-j - Rᵢ = _sylpoly_shift(s, i+1, u) - Rⱼ = _sylpoly(s, j+1, v) - Rₖ = _sylpoly_shift(s, k+1, w) + k = (Degree+2)-i-j + Rᵢ = gen_sylpoly_shift(s, i+1, u, T) + Rⱼ = gen_sylpoly(s, j+1, v, T) + Rₖ = gen_sylpoly_shift(s, k+1, w, T) - 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ₖ - - nts[idx] = (value=Rᵢ*Rⱼ*Rₖ*nd2, curl=curl) - idx += 1 + dRᵢ = gen_sylpoly_shift_diff(s, i+1, u, T) + dRⱼ = gen_sylpoly_diff(s, j+1, v, T) + dRₖ = gen_sylpoly_shift_diff(s, k+1, w, T) + 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ₖ ) + + nts = :($nts..., (value=$Rᵢ*$Rⱼ*$Rₖ*$nd2, curl=$curl)) end - for i in 1:d+1 - j = (d+2)-i + for i in 1:Degree+1 + j = (Degree+2)-i k = 0 - Rᵢ = _sylpoly_shift(s, i+1, u) - Rⱼ = _sylpoly_shift(s, j+1, v) - Rₖ = _sylpoly(s, k+1, w) + Rᵢ = gen_sylpoly_shift(s, i+1, u, T) + Rⱼ = gen_sylpoly_shift(s, j+1, v, T) + Rₖ = gen_sylpoly(s, k+1, w, T) - dRᵢ = _sylpoly_shift_diff(s, i+1, u) - dRⱼ = _sylpoly_shift_diff(s, j+1, v) - dRₖ = _sylpoly_diff(s, k+1, w) + dRᵢ = gen_sylpoly_shift_diff(s, i+1, u, T) + dRⱼ = gen_sylpoly_shift_diff(s, j+1, v, T) + dRₖ = gen_sylpoly_diff(s, k+1, w, T) - 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ₖ + 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ₖ ) - nts[idx] = (value=Rᵢ*Rⱼ*Rₖ*nd3, curl=curl) - idx += 1 + nts = :($nts..., (value=$Rᵢ*$Rⱼ*$Rₖ*$nd3, curl=$curl)) end - for i in 1:d+1 - for j in 1:d+1 - k = (d+2)-i-j + for i in 1:Degree+1 + for j in 1:Degree+1 + k = (Degree+2)-i-j 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) - S1 = Rᵢ*Rsⱼ*Rsₖ*nd1 - S2 = Rsᵢ*Rⱼ*Rsₖ*nd2 - S3 = Rsᵢ*Rsⱼ*Rₖ*nd3 - - 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ₖ - - nts[idx] = (value=S2-S3, curl=curlS2 - curlS3) - idx += 1 + Rsᵢ = gen_sylpoly_shift(s, i+1, u, T) + Rsⱼ = gen_sylpoly_shift(s, j+1, v, T) + Rsₖ = gen_sylpoly_shift(s, k+1, w, T) + Rᵢ = gen_sylpoly(s, i+1, u, T) + Rⱼ = gen_sylpoly(s, j+1, v, T) + Rₖ = gen_sylpoly(s, k+1, w, T) + S1 = :( $Rᵢ*$Rsⱼ*$Rsₖ*$nd1 ) + S2 = :( $Rsᵢ*$Rⱼ*$Rsₖ*$nd2 ) + S3 = :( $Rsᵢ*$Rsⱼ*$Rₖ*$nd3 ) + + dRsᵢ = gen_sylpoly_shift_diff(s, i+1, u, T) + dRsⱼ = gen_sylpoly_shift_diff(s, j+1, v, T) + dRsₖ = gen_sylpoly_shift_diff(s, k+1, w, T) + dRᵢ = gen_sylpoly_diff(s, i+1, u, T) + dRⱼ = gen_sylpoly_diff(s, j+1, v, T) + dRₖ = gen_sylpoly_diff(s, k+1, w, T) + + 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ₖ ) + + nts = :($nts..., (value=$S2-$S3, curl=$curlS2 - $curlS3)) - nts[idx] = (value=S3-S1, curl=curlS3 - curlS1) - idx += 1 + nts = :($nts..., (value=$S3-$S1, curl=$curlS3 - $curlS1)) end end - return SVector{NF}(nts) + ex = :(SVector{NF}(())) + + for i in 1:NF + push!(ex.args[2].args, :($nts[$i])) + end + + return ex end diff --git a/src/utils/lagpolys.jl b/src/utils/lagpolys.jl index 60620706..e01fa77b 100644 --- a/src/utils/lagpolys.jl +++ b/src/utils/lagpolys.jl @@ -30,18 +30,45 @@ function _lagpoly_diff(nodes, i, s, i0=1, i1=length(nodes)) return r end -function _sylpoly(nodes, i, s) - _lagpoly(nodes, i, s, 1, i) + +# Versions of the above functions to be used with @generated functions. These return expressions +# that can be evaluated at compile time, which allows for more efficient code when the number of +# nodes is known at compile time. +function gen_lagpoly(nodes, i, s, i0, i1, T) + p = :(one(T)) + for j in i0:i1 + j == i && continue + p = :(($p * ($s - $nodes[$j]) / ($nodes[$i] - $nodes[$j]))) + end + return p +end + +function gen_lagpoly_diff(nodes, i, s, i0, i1, T) + dp = :(zero(T)) + for j in i0:i1 + j == i && continue + p = :(one(T)) + for k in i0:i1 + (k == i || k == j) && continue + p = :(($p * ($s - $nodes[$k]) / ($nodes[$i] - $nodes[$k]))) + end + dp = :(($dp + $p / ($nodes[$i] - $nodes[$j]))) + end + return dp +end + +function gen_sylpoly(nodes, i, s, T) + gen_lagpoly(nodes, i, s, 1, i, T) end -function _sylpoly_diff(nodes, i, s) - _lagpoly_diff(nodes, i, s, 1, i) +function gen_sylpoly_diff(nodes, i, s, T) + gen_lagpoly_diff(nodes, i, s, 1, i, T) end -function _sylpoly_shift(nodes, i, s) - _lagpoly(nodes, i, s, 2, i) +function gen_sylpoly_shift(nodes, i, s, T) + gen_lagpoly(nodes, i, s, 2, i, T) end -function _sylpoly_shift_diff(nodes, i, s) - _lagpoly_diff(nodes, i, s, 2, i) +function gen_sylpoly_shift_diff(nodes, i, s, T) + gen_lagpoly_diff(nodes, i, s, 2, i, T) end \ No newline at end of file From 039fcac96ce232417b36fc5068d76a271b9c51d9 Mon Sep 17 00:00:00 2001 From: Cedric Muenger Date: Mon, 4 May 2026 18:29:18 +0200 Subject: [PATCH 2/6] generated Lagrange shapefunctions --- src/BEAST.jl | 3 +- src/bases/global/globaldofs.jl | 4 +- src/bases/lagrange.jl | 2036 +++++++++++------- src/bases/local/gwpdivlocal.jl | 2 +- src/bases/local/laglocal.jl | 738 ++++--- test/test_basis.jl | 16 +- test/test_gridfunction.jl | 10 +- test/test_higher_order_lagrange_functions.jl | 6 +- test/test_local_assembly.jl | 2 +- test/test_trace.jl | 2 +- 10 files changed, 1743 insertions(+), 1076 deletions(-) diff --git a/src/BEAST.jl b/src/BEAST.jl index bde48f83..8b8d9861 100644 --- a/src/BEAST.jl +++ b/src/BEAST.jl @@ -160,6 +160,8 @@ include("utils/lagpolys.jl") include("utils/butchertableau.jl") include("utils/variational.jl") +include("bases/global/globaldofs.jl") + include("bases/localbasis.jl") include("bases/local/interpolate.jl") @@ -183,7 +185,6 @@ include("bases/lincomb.jl") 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") diff --git a/src/bases/global/globaldofs.jl b/src/bases/global/globaldofs.jl index ad7908bd..e58da5e7 100644 --- a/src/bases/global/globaldofs.jl +++ b/src/bases/global/globaldofs.jl @@ -126,14 +126,14 @@ end @test size(dofs) == (10,2) A = [ 0 0 - 0 1 - 1 0 0 0 0 0 0 0 0 0 0 0 + 1 0 0 0 + 0 1 0 0] @test dofs ≈ A end \ No newline at end of file diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl index ef26ac56..36aa7b78 100644 --- a/src/bases/lagrange.jl +++ b/src/bases/lagrange.jl @@ -29,26 +29,558 @@ end refspace(space::LagrangeBasis{D,C,M,T,NF}) where {D,C,M,T,NF} = LagrangeRefSpace{T,D,dimension(geometry(space))+1,NF}() subset(s::S,I) where {S <: Space} = S(s.geo, s.fns[I], s.pos[I]) -function lagrangecxd0(mesh) - U = universedimension(mesh) - D1 = dimension(mesh)+1 +#==== +LAGRANGE SPACES +====# + +lagrangecxd0(mesh) = lagrangecx(mesh; order=0) +lagrangec0d1(mesh; dirichlet=false) = lagrangec0(mesh; dirichlet=dirichlet, order=1) +lagrangec0d2(mesh; dirichlet=false) = lagrangec0(mesh; dirichlet=dirichlet, order=2) + +lagrangec0d1(mesh, bnd) = lagrangec0(mesh, bnd; order=1) + +#==== +LAGRANGE SPACES on SEGMENTS +====# +function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{<:Any, 2}; dirichlet::Bool=false, order=1) + + @assert order>0 + if dirichlet == false + # return lagrangec0d1(mesh, boundary(mesh)) + return lagrangec0(mesh, skeleton(mesh,0); order=order) + else + + verts = skeleton(mesh, 0) + bnd = boundary(mesh) + bndverts = skeleton(bnd, 0) + + overlap = overlap_gpredicate(bndverts) + notoverlap = (m,c) -> !overlap(chart(m,c)) + interiorverts = submesh(notoverlap,verts) + + return lagrangec0(mesh, interiorverts; order=order) + end +end + +# build continuous Lagrange elements on a 1D mesh +function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{<:Any,2}, + nodes::CompScienceMeshes.AbstractMesh{<:Any,1}; order=1) + + if order > 10 + error("Maximal 1D Lagrange polynomial order supported is 10") + end + + Conn = connectivity(nodes, mesh, abs) + rows = rowvals(Conn) + vals = nonzeros(Conn) + T = coordtype(mesh) - geometry = mesh - num_cells = numcells(mesh) + P = vertextype(mesh) + S = Shape{T} - # create the local shapes - fns = Vector{Vector{Shape{T}}}(undef,num_cells) - pos = Vector{vertextype(mesh)}(undef,num_cells) - for (i,cell) in enumerate(mesh) - fns[i] = [Shape(i, 1, T(1.0))] - pos[i] = cartesian(center(chart(mesh, cell))) - end + fns = Vector{Vector{S}}() + pos = Vector{P}() + + u = T(1.0) + + for (i,node) in enumerate(nodes) + fn = Vector{S}() + for k in nzrange(Conn,i) + cellid = rows[k] + refid = vals[k] + push!(fn, Shape(cellid, refid, u)) + end + push!(fns,fn) + push!(pos,cartesian(center(chart(nodes,node)))) + end + + NF = order + 1 + + for (c, cell) in enumerate(mesh) + ch = chart(mesh, cell) + for r in 3:NF + push!(fns, S[S(c, r, u)]) + push!(pos, cartesian(neighborhood(ch, SVector(u - (r-2)/(NF-1))))) + end + end + + LagrangeBasis{order,0,NF}(mesh, fns, pos) +end + +# build discontinuous Lagrange elements on a 1D mesh +function lagrangecx(mesh::CompScienceMeshes.AbstractMesh{<:Any,2}; order) + + T = coordtype(mesh) + + P = vertextype(mesh) + S = Shape{T} + + fns = Vector{Vector{S}}() + pos = Vector{P}() + + u = one(T) + + if order == 0 + for (c,cell) in enumerate(mesh) + ch = chart(mesh,cell) + + push!(fns, S[S(c,1,u)]) + push!(pos, cartesian(neighborhood(ch, T(0.5)))) + end + else + for (c,cell) in enumerate(mesh) + ch = chart(mesh,cell) + + push!(fns, S[S(c,1,u)]) + push!(pos, cartesian(neighborhood(ch, T(1.0)))) + + push!(fns, S[S(c,2,u)]) + push!(pos, cartesian(neighborhood(ch, T(0.0)))) + end + end + + NF = order + 1 + + for (c,cell) in enumerate(mesh) + ch = chart(mesh,cell) + for r in 3:NF + push!(fns, S[S(c,r,u)]) + push!(pos, cartesian(neighborhood(ch, u - (r-2)/(NF-1)))) + end + end + + return LagrangeBasis{order,-1,NF}(mesh, fns, pos) +end + + +#==== +LAGRANGE SPACES on TRIANGLES +====# +function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{<:Any, 3}; dirichlet::Bool=false, order=1) + + @assert order>0 + if dirichlet == false + # return lagrangec0d1(mesh, boundary(mesh)) + return lagrangec0(mesh, skeleton(mesh,0), skeleton(mesh,1); order=order) + else + + bnd = boundary(mesh) + + verts = skeleton(mesh, 0) + bndverts = skeleton(bnd, 0) + overlap = overlap_gpredicate(bndverts) + notoverlap = (m,c) -> !overlap(chart(m,c)) + interiorverts = submesh(notoverlap,verts) + + edges = skeleton(mesh, 1) + bndeges = skeleton(bnd, 1) + overlap = overlap_gpredicate(bndeges) + notoverlap = (m,c) -> !overlap(chart(m,c)) + interioredges = submesh(notoverlap,edges) + + return lagrangec0(mesh, interiorverts, interioredges; order=order) + end +end + +function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{U,3}, + neumannbnd::CompScienceMeshes.AbstractMesh{U,2}; order=1) where {U} + + bnd = boundary(mesh) + + # Remove vertices at the boundary of the neumann boundary + bndverts = skeleton(boundary(neumannbnd), 0) + neumannbndverts = skeleton(neumannbnd, 0) + overlap = overlap_gpredicate(bndverts) + notoverlap = (m,c) -> !overlap(chart(m,c)) + interiorneumannbndverts = submesh(notoverlap,neumannbndverts) + + # Remove the interior neumann vertices from the boundary vertices + # to obtain the dirichlet vertices on the boundary + bndverts = skeleton(bnd, 0) + overlap = overlap_gpredicate(interiorneumannbndverts) + notoverlap = (m,c) -> !overlap(chart(m,c)) + dirichletbndverts = submesh(notoverlap,bndverts) + + @show dirichletbndverts + # Remove dirchlet vertices from the mesh vertices + # to obtain the interior vertices + verts = skeleton(mesh, 0) + overlap = overlap_gpredicate(dirichletbndverts) + vertoverlap = (m,c) -> !overlap(chart(m,c)) + interiorverts = submesh(vertoverlap,verts) + + # Remove neumann edges from the boundary edges + # to obtain the dirichlet edges + neumannbndedges = skeleton(neumannbnd, 1) + bndedges = skeleton(bnd, 1) + overlap = overlap_gpredicate(neumannbndedges) + notoverlap = (m,c) -> !overlap(chart(m,c)) + dirichletbndedges = submesh(notoverlap,bndedges) + + # Remove dirichlet edges from the mesh edges + # to obtain the interior edges + edges = skeleton(mesh, 1) + overlap = overlap_gpredicate(dirichletbndedges) + edgeoverlap = (m,c) -> !overlap(chart(m,c)) + interioredges = submesh(edgeoverlap,edges) + + return lagrangec0(mesh, interiorverts, interioredges; order=order) +end + + +# build continuous Lagrange elements on a 2D mesh +function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{U,3}, + nodes::CompScienceMeshes.AbstractMesh{U,1}, + edges::CompScienceMeshes.AbstractMesh{U,2}; order=1) where {U} + + T = coordtype(mesh) + atol = sqrt(eps(T)) + + P = vertextype(mesh) + S = Shape{T} + F = Vector{S} + + # nodes = skeleton(mesh, 0) + # edges = skeleton(mesh, 1) + + C02 = connectivity(mesh, nodes, abs); R02 = rowvals(C02); V02 = nonzeros(C02) + C12 = connectivity(mesh, edges, abs); R12 = rowvals(C12); V12 = nonzeros(C12) + + ne = order-1 + nf = div((order-2)*(order-1), 2) + + nV = length(nodes) + nE = length(edges) * ne + nF = length(mesh) * nf + N = nV + nE + nF + + localspace = LagrangeRefSpace{T,order,3,binomial(2+order,2)}() + dom = domain(chart(mesh, first(mesh))) + localdim = numfunctions(localspace, dom) + + d = order + fns = [S[] for n in 1:(nV+nE+nF)] + pos = fill(point(0,0,0), nV+nE+nF) + for cell in mesh + cell_ch = chart(mesh, cell) + V = R02[nzrange(C02,cell)] + I = V02[nzrange(C02,cell)] + for (i,v) in zip(I, V) + vertex_ch = chart(nodes, v) + gids = [v] + lids = localindices(_LagrangeGlobalNodesDoFs(d), cell_ch, localspace, i) + v = globaldofs(vertex_ch, cell_ch, localspace, _LagrangeGlobalNodesDoFs(d)) + α = v[lids, :] + β = inv(α') + for i in axes(β,1) + for j in axes(β,2) + isapprox(β[i,j], 0; atol) && continue + push!(fns[gids[i]], S(cell, lids[j], β[i,j])) + end + end + end + + order < 2 && continue + E = R12[nzrange(C12,cell)] + I = V12[nzrange(C12,cell)] + for (i,e) in zip(I, E) + edge_ch = chart(edges, e) + gids = nV + (e-1)*ne + 1: nV + e*ne + lids = localindices(_LagrangeGlobalEdgeDoFs(d), cell_ch, localspace, i) + @assert length(lids) == length(gids) + v = globaldofs(edge_ch, cell_ch, localspace, _LagrangeGlobalEdgeDoFs(d)) + α = v[lids, :] + β = inv(α') + for i in axes(β,1) + for j in axes(β,2) + isapprox(β[i,j], 0; atol) && continue + push!(fns[gids[i]], S(cell, lids[j], β[i,j])) + end + end + end + + order < 3 && continue + F = [cell] + for (q,f) in enumerate(F) + face_ch = chart(mesh, f) + gids = nV + nE + (f-1)*nf + 1 : nV + nE + f*nf + lids = localindices(_LagrangeGlobalFaceDoFs(d), cell_ch, localspace, 1) + v = globaldofs(face_ch, cell_ch, localspace, _LagrangeGlobalFaceDoFs(d)) + α = v[lids, :] + # @show α + β = inv(α') + for i in axes(β,1) + for j in axes(β,2) + isapprox(β[i,j], 0; atol) && continue + push!(fns[gids[i]], S(cell, lids[j], β[i,j])) + end + end + end + end + + for v in nodes pos[v] = cartesian(center(chart(nodes, v))) end + for e in edges pos[nV + (e-1)*ne + 1: nV + e*ne] .= Ref(cartesian(center(chart(edges, e)))) end + for f in mesh pos[nV + nE + (f-1)*nf + 1: nV + nE + f*nf] .= Ref(cartesian(center(chart(mesh, f)))) end - NF = 1 - LagrangeBasis{0,-1,NF}(geometry, fns, pos) + return LagrangeBasis{order,0,localdim}(mesh, fns, pos) +end + + +# build discontinuous Lagrange elements on a 2D mesh +function lagrangecx(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}; order) + + T = coordtype(mesh) + NF = binomial(2+order, 2) + P = vertextype(mesh) + S = Shape{T} + + fns = Vector{Vector{S}}(undef, length(mesh) * NF) + pos = Vector{P}(undef, length(mesh) * NF) + + idx = 1 + u = one(T) + for (c,cell) in enumerate(mesh) + ch = chart(mesh,cell) + for r in 1:NF + fns[idx] = S[S(c,r,u)] + #TODO: use actual node position instead of center of element + pos[idx] = cartesian(center(ch)) + idx += 1 + end + end + + return LagrangeBasis{order,-1,NF}(mesh, fns, pos) +end + + + +#==== +LAGRANGE SPACES on TETRAHEDRA +====# +function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{<:Any, 4}; dirichlet::Bool=true, order=1) + + @assert order>0 + if dirichlet == false + # return lagrangec0d1(mesh, boundary(mesh)) + return lagrangec0(mesh, skeleton(mesh,0), skeleton(mesh,1), skeleton(mesh,2); order=order) + else + + bnd = boundary(mesh) + + verts = skeleton(mesh, 0) + bndverts = skeleton(bnd, 0) + overlap = overlap_gpredicate(bndverts) + notoverlap = (c,i) -> !overlap(c,i) + interiorverts = submesh(notoverlap,verts) + + edges = skeleton(mesh, 1) + bndeges = skeleton(bnd, 1) + overlap = overlap_gpredicate(bndeges) + notoverlap = (c,i) -> !overlap(c,i) + interioredges = submesh(notoverlap,edges) + + + faces = skeleton(mesh, 2) + bndfaces = skeleton(bnd, 2) + overlap = overlap_gpredicate(bndfaces) + notoverlap = (c,i) -> !overlap(c,i) + interiorfaces = submesh(notoverlap,faces) + + return lagrangec0(mesh, interiorverts, interioredges, interiorfaces; order=order) + end +end + + + + +# build continuous Lagrange elements on a 3D mesh +function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{U,4}, + nodes::CompScienceMeshes.AbstractMesh{U,1}, + edges::CompScienceMeshes.AbstractMesh{U,2}, + faces::CompScienceMeshes.AbstractMesh{U,3}; order=1) where {U} + + #TODO generalize for higher order + if order > 1 + error("Only linear 3D Lagrange elements are supported") + end + + T = coordtype(mesh) + atol = sqrt(eps(T)) + + S = Shape{T} + + C02 = connectivity(mesh, nodes, abs); R02 = rowvals(C02); V02 = nonzeros(C02) + + nV = length(nodes) + + localspace = LagrangeRefSpace{T,order,4,4}() + dom = domain(chart(mesh, first(mesh))) + localdim = numfunctions(localspace, dom) + + d = order + fns = [S[] for n in 1:(nV)] + pos = fill(point(0,0,0), nV) + for cell in mesh + cell_ch = chart(mesh, cell) + V = R02[nzrange(C02,cell)] + I = V02[nzrange(C02,cell)] + for (i,v) in zip(I, V) + vertex_ch = chart(nodes, v) + gids = [v] + lids = localindices(_LagrangeGlobalNodesDoFs(d), cell_ch, localspace, i) + v = globaldofs(vertex_ch, cell_ch, localspace, _LagrangeGlobalNodesDoFs(d)) + α = v[lids, :] + β = inv(α') + for i in axes(β,1) + for j in axes(β,2) + isapprox(β[i,j], 0; atol) && continue + push!(fns[gids[i]], S(cell, lids[j], β[i,j])) + end + end + end + + + end + + for v in nodes pos[v] = cartesian(center(chart(nodes, v))) end + + return LagrangeBasis{order,0,localdim}(mesh, fns, pos) +end + + +# build discontinuous Lagrange elements on a 3D mesh +function lagrangecx(mesh::CompScienceMeshes.AbstractMesh{<:Any,4}; order) + + #TODO generalize for higher order + if order > 1 + error("Only up to linear 3D Lagrange elements are supported") + end + + T = coordtype(mesh) + NF = binomial(2+order, 2) + P = vertextype(mesh) + S = Shape{T} + + fns = Vector{Vector{S}}(undef, length(mesh) * NF) + pos = Vector{P}(undef, length(mesh) * NF) + + idx = 1 + u = one(T) + for (c,cell) in enumerate(mesh) + ch = chart(mesh,cell) + for r in 1:NF + fns[idx] = S[S(c,r,u)] + #TODO: use actual node position instead of center of element + pos[idx] = cartesian(center(ch)) + idx += 1 + end + end + + return LagrangeBasis{order,-1,NF}(mesh, fns, pos) end +gradient(space::LagrangeBasis{1,0,<:CompScienceMeshes.AbstractMesh{<:Any,2}}, geo, fns) = + LagrangeBasis{0,-1,1}(geo, fns, space.pos) + +gradient(space::LagrangeBasis{1,0,<:CompScienceMeshes.AbstractMesh{<:Any,3}}, geo, fns) = + NDBasis(geo, fns, space.pos) + +gradient(space::LagrangeBasis{1,0,<:CompScienceMeshes.AbstractMesh{<:Any,4}}, geo, fns) = + NDLCCBasis(geo, fns) + + +function curl(X::LagrangeBasis{D,C,M,T,NF,P} , geo, fns) where {D,C,M,T,NF,P} + GWPDivSpace{T,M,Vector{P}}(geo, fns, deepcopy(positions(X)),D-1) +end + + +@testitem "curl - global" begin + using LinearAlgebra + using CompScienceMeshes + const CSM = CompScienceMeshes + + T = Float64 + + m = CSM.meshrectangle(1.0, 1.0, 0.5, 3) + X = BEAST.lagrangec0(m; order=1) + curlX = BEAST.curl(X) + + x = BEAST.refspace(X) + curlx = BEAST.refspace(curlX) + + 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 = zeros(T,3) + ϕp = x(p) + for sh in X.fns[i] + sh.cellid == cellid || continue + r1 += sh.coeff * ϕp[sh.refid].curl + end + + r2 = zeros(T,3) + ϕp = curlx(p) + for sh in curlX.fns[i] + sh.cellid == cellid || continue + r2 += sh.coeff * ϕp[sh.refid].value + end + global err = max(err, norm(r1-r2)) + end + end + + @test err < sqrt(eps(T)) +end + +# +# Sclar trace for Laggrange element based spaces +# +function strace(X::LagrangeBasis{1,0}, geo, fns::Vector) + # dimension of the boundary + n = dimension(geo) + # degree of the space + d = 1 + # number of local shape functions + NF = binomial(n+d,n) + + trpos = deepcopy(positions(X)) + LagrangeBasis{1,0,NF}(geo, fns, trpos) +end + + + +# function lagrangecxd0(mesh) + +# U = universedimension(mesh) +# D1 = dimension(mesh)+1 +# T = coordtype(mesh) +# geometry = mesh +# num_cells = numcells(mesh) + +# # create the local shapes +# fns = Vector{Vector{Shape{T}}}(undef,num_cells) +# pos = Vector{vertextype(mesh)}(undef,num_cells) +# for (i,cell) in enumerate(mesh) +# fns[i] = [Shape(i, 1, T(1.0))] +# pos[i] = cartesian(center(chart(mesh, cell))) +# end + +# NF = 1 +# LagrangeBasis{0,-1,NF}(geometry, fns, pos) +# end + + +#==== +UNIT FUNCTION +====# + """ unitfunctioncxd0(mesh) @@ -176,51 +708,65 @@ function unitfunctionc0d1(mesh, nodes::CompScienceMeshes.AbstractMesh{U,1} where LagrangeBasis{1,0,NF}(mesh, fns, pos) end -""" - lagrangec0d1(mesh[, bnd]) -Construct the basis of continuous, piecewise linear basis functions subordinate -to mesh `mesh`. Basis functions are constructed at vertices in the interior of -the mesh and on the closure of 'bnd'. In particular, leaving out the second argument -creates a finite element space subject to homogeneous Dirichlet boundary conditions. -""" -function lagrangec0d1_dirichlet(mesh) +function _surface(X) + C = unitfunctioncxd0(X.geo) + G = assemble(Identity(),X,X) + u = Vector(assemble(Identity(),X,C)[1:end,1]) + dot(u,G\u) +end + +# """ +# lagrangec0d1(mesh[, bnd]) + +# Construct the basis of continuous, piecewise linear basis functions subordinate +# to mesh `mesh`. Basis functions are constructed at vertices in the interior of +# the mesh and on the closure of 'bnd'. In particular, leaving out the second argument +# creates a finite element space subject to homogeneous Dirichlet boundary conditions. +# """ +# function lagrangec0d1_dirichlet(mesh) - verts = skeleton(mesh, 0) - detached = trues(numvertices(mesh)) - for v in cells(verts) - detached[v.indices] = false - end +# verts = skeleton(mesh, 0) +# detached = trues(numvertices(mesh)) +# for v in cells(verts) +# detached[v.indices] = false +# end - bnd = boundary(mesh) - bndverts = skeleton(bnd, 0) - notonbnd = trues(numvertices(mesh)) - for v in cells(bndverts) - notonbnd[v.indices] = false - end +# bnd = boundary(mesh) +# bndverts = skeleton(bnd, 0) +# notonbnd = trues(numvertices(mesh)) +# for v in cells(bndverts) +# notonbnd[v.indices] = false +# end - vertexlist = findall(notonbnd .& .!detached) - lagrangec0d1(mesh, vertexlist, Val{dimension(mesh)+1}) -end +# vertexlist = findall(notonbnd .& .!detached) +# lagrangec0d1(mesh, vertexlist, Val{dimension(mesh)+1}) +# end -function lagrangec0_dirichlet(mesh; order=1) +# function lagrangec0_dirichlet(mesh; order=1) - verts = skeleton(mesh, 0) - detached = trues(numvertices(mesh)) - for v in cells(verts) - detached[v.indices[1]] = false - end +# verts = skeleton(mesh, 0) +# detached = trues(numvertices(mesh)) +# for v in cells(verts) +# detached[v.indices[1]] = false +# end - bnd = boundary(mesh) - bndverts = skeleton(bnd, 0) - notonbnd = trues(numvertices(mesh)) - for v in cells(bndverts) - notonbnd[v.indices[1]] = false - end +# bnd = boundary(mesh) +# bndverts = skeleton(bnd, 0) +# notonbnd = trues(numvertices(mesh)) +# for v in cells(bndverts) +# notonbnd[v.indices[1]] = false +# end + +# vertexlist = findall(notonbnd .& .!detached) +# lagrangec0(mesh, vertexlist, Val{dimension(mesh)+1}; order=order) +# end + + +#==== +DUAL LAGRANGE SPACES +====# - vertexlist = findall(notonbnd .& .!detached) - lagrangec0(mesh, vertexlist, Val{dimension(mesh)+1}; order=order) -end function interior_and_junction_vertices(mesh, jct) verts = skeleton(mesh, 0) @@ -331,296 +877,296 @@ function singleduallagd0(fine, F, v; interpolatory=false) return fn end -""" - lagrangec0d1(mesh; dirichlet=[true|false]) -> basis - -Build lagrangec0d1 elements, including (dirichlet=false) -or excluding (dirichlet=true) those attached to boundary vertices. -""" -function lagrangec0d1(mesh; dirichlet::Bool=true) - if dirichlet == false - # return lagrangec0d1(mesh, boundary(mesh)) - return lagrangec0d1(mesh, skeleton(mesh,0)) - else - return lagrangec0d1_dirichlet(mesh) - end -end - -function lagrangec0(mesh; order=1, dirichlet::Bool=true) - - @assert order>0 - if dirichlet == false - # return lagrangec0d1(mesh, boundary(mesh)) - return lagrangec0(mesh, skeleton(mesh,0); order=order) - else - return lagrangec0_dirichlet(mesh; order=order) - end -end - -function lagrangec0d1(mesh, jct) - vertexlist = interior_and_junction_vertices(mesh, jct) - lagrangec0d1(mesh, vertexlist, Val{dimension(mesh)+1}) -end - -# build continuous linear Lagrange elements on a 2D manifold -function lagrangec0d1(mesh, vertexlist::Vector, ::Type{Val{3}}) - - T = coordtype(mesh) - U = universedimension(mesh) - - cellids, ncells = vertextocellmap(mesh) - - Cells = cells(mesh) - Verts = vertices(mesh) - - # create the local shapes - fns = Vector{Shape{T}}[] - pos = Vector{vertextype(mesh)}() - - sizehint!(fns, length(vertexlist)) - sizehint!(pos, length(vertexlist)) - for v in vertexlist - - numshapes = ncells[v] - numshapes == 0 && continue - - shapes = Vector{Shape{T}}(undef,numshapes) - for s in 1: numshapes - c = cellids[v,s] - # cell = mesh.faces[c] - cell = Cells[c] - - localid = something(findfirst(isequal(v), cell.indices),0) - @assert localid != 0 - - shapes[s] = Shape(c, localid, T(1.0)) - end - - push!(fns, shapes) - push!(pos, Verts[v]) - end - - NF = 3 - LagrangeBasis{1,0,NF}(mesh, fns, pos) -end - -# for manifolds of dimension 1 -function lagrangec0d1(mesh, vertexlist, ::Type{Val{2}}) - - T = coordtype(mesh) - U = universedimension(mesh) - P = vertextype(mesh) - - geometry = mesh - - cellids, ncells = vertextocellmap(mesh) - - # create the local shapes - numverts = numvertices(mesh) - - fns = Vector{Vector{Shape{T}}}() - pos = Vector{P}() - - sizehint!(fns, length(vertexlist)) - sizehint!(pos, length(vertexlist)) - for v in vertexlist - - numshapes = ncells[v] - numshapes == 0 && continue # skip detached vertices - - shapes = Vector{Shape{T}}(undef,numshapes) - for s in 1: numshapes - c = cellids[v,s] - cell = mesh.faces[c] - if cell[1] == v - shapes[s] = Shape(c, 1, T(1.0)) - elseif cell[2] == v - shapes[s] = Shape(c, 2, T(1.0)) - else - error("Junctions not supported") - end - end - - push!(fns, shapes) - push!(pos, mesh.vertices[v]) - end - - NF = 2 - LagrangeBasis{1,0,NF}(geometry, fns, pos) -end - -function lagrangec0(mesh, vertexlist, ::Type{Val{2}}; order=1) - - if order > 10 - error("Maximal 1D Lagrange polynomial order supported is 10") - end - - T = coordtype(mesh) - U = universedimension(mesh) - P = vertextype(mesh) - S = Shape{T} - - geometry = mesh - - cellids, ncells = vertextocellmap(mesh) - - # create the local shapes - numverts = numvertices(mesh) - - fns = Vector{Vector{S}}() - pos = Vector{P}() - - sizehint!(fns, length(vertexlist)) - sizehint!(pos, length(vertexlist)) - for v in vertexlist - - numshapes = ncells[v] - numshapes == 0 && continue # skip detached vertices - - shapes = Vector{S}(undef,numshapes) - for s in 1: numshapes - c = cellids[v,s] - cell = mesh.faces[c] - if cell[1] == v - shapes[s] = Shape(c, 1, T(1.0)) - elseif cell[2] == v - shapes[s] = Shape(c, 2, T(1.0)) - else - error("Junctions not supported") - end - end - - push!(fns, shapes) - push!(pos, mesh.vertices[v]) - end - - NF = order + 1 - - for (c,cell) in enumerate(mesh) - ch = chart(mesh,cell) - for r in 3:NF - push!(fns, S[S(c,r,T(1.0))]) - push!(pos, cartesian(neighborhood(ch, (r-2)/(NF-1)))) - end - end - - LagrangeBasis{order,0,NF}(geometry, fns, pos) -end - -function lagrangec0d1(mesh, nodes::CompScienceMeshes.AbstractMesh{<:Any,1}) - - Conn = connectivity(nodes, mesh, abs) - rows = rowvals(Conn) - vals = nonzeros(Conn) - - T = coordtype(mesh) - P = vertextype(mesh) - S = Shape{T} - - fns = Vector{Vector{S}}() - pos = Vector{P}() - for (i,node) in enumerate(nodes) - fn = Vector{S}() - for k in nzrange(Conn,i) - cellid = rows[k] - refid = vals[k] - push!(fn, Shape(cellid, refid, T(1.0))) - end - push!(fns,fn) - push!(pos,cartesian(center(chart(nodes,node)))) - end - - NF = dimension(mesh) + 1 - - LagrangeBasis{1,0,NF}(mesh, fns, pos) -end +# """ +# lagrangec0d1(mesh; dirichlet=[true|false]) -> basis -function lagrangec0(mesh, nodes::CompScienceMeshes.AbstractMesh{<:Any,1}; order=1) +# Build lagrangec0d1 elements, including (dirichlet=false) +# or excluding (dirichlet=true) those attached to boundary vertices. +# """ +# function lagrangec0d1(mesh; dirichlet::Bool=true) +# if dirichlet == false +# # return lagrangec0d1(mesh, boundary(mesh)) +# return lagrangec0d1(mesh, skeleton(mesh,0)) +# else +# return lagrangec0d1_dirichlet(mesh) +# end +# end + +# function lagrangec0(mesh; order=1, dirichlet::Bool=true) + +# @assert order>0 +# if dirichlet == false +# # return lagrangec0d1(mesh, boundary(mesh)) +# return lagrangec0(mesh, skeleton(mesh,0); order=order) +# else +# return lagrangec0_dirichlet(mesh; order=order) +# end +# end - if order > 10 - error("Maximal 1D Lagrange polynomial order supported is 10") - end +# function lagrangec0d1(mesh, jct) +# vertexlist = interior_and_junction_vertices(mesh, jct) +# lagrangec0d1(mesh, vertexlist, Val{dimension(mesh)+1}) +# end - Conn = connectivity(nodes, mesh, abs) - rows = rowvals(Conn) - vals = nonzeros(Conn) +# # build continuous linear Lagrange elements on a 2D manifold +# function lagrangec0d1(mesh, vertexlist::Vector, ::Type{Val{3}}) - T = coordtype(mesh) - P = vertextype(mesh) - S = Shape{T} +# T = coordtype(mesh) +# U = universedimension(mesh) - fns = Vector{Vector{S}}() - pos = Vector{P}() +# cellids, ncells = vertextocellmap(mesh) - u = T(1.0) +# Cells = cells(mesh) +# Verts = vertices(mesh) - for (i,node) in enumerate(nodes) - fn = Vector{S}() - for k in nzrange(Conn,i) - cellid = rows[k] - refid = vals[k] - push!(fn, Shape(cellid, refid, u)) - end - push!(fns,fn) - push!(pos,cartesian(center(chart(nodes,node)))) - end +# # create the local shapes +# fns = Vector{Shape{T}}[] +# pos = Vector{vertextype(mesh)}() - NF = order + 1 +# sizehint!(fns, length(vertexlist)) +# sizehint!(pos, length(vertexlist)) +# for v in vertexlist - for (c, cell) in enumerate(mesh) - ch = chart(mesh, cell) - for r in 3:NF - push!(fns, S[S(c, r, u)]) - push!(pos, cartesian(neighborhood(ch, SVector(u - (r-2)/(NF-1))))) - end - end +# numshapes = ncells[v] +# numshapes == 0 && continue - LagrangeBasis{order,0,NF}(mesh, fns, pos) -end +# shapes = Vector{Shape{T}}(undef,numshapes) +# for s in 1: numshapes +# c = cellids[v,s] +# # cell = mesh.faces[c] +# cell = Cells[c] + +# localid = something(findfirst(isequal(v), cell.indices),0) +# @assert localid != 0 + +# shapes[s] = Shape(c, localid, T(1.0)) +# end -function lagrangec0d2(mesh::CompScienceMeshes.AbstractMesh{U,3}, - nodes::CompScienceMeshes.AbstractMesh{U,1}, - edges::CompScienceMeshes.AbstractMesh{U,2}) where {U} +# push!(fns, shapes) +# push!(pos, Verts[v]) +# end + +# NF = 3 +# LagrangeBasis{1,0,NF}(mesh, fns, pos) +# end + +# # for manifolds of dimension 1 +# function lagrangec0d1(mesh, vertexlist, ::Type{Val{2}}) + +# T = coordtype(mesh) +# U = universedimension(mesh) +# P = vertextype(mesh) + +# geometry = mesh - Conn = connectivity(nodes, mesh, abs) - rows = rowvals(Conn) - vals = nonzeros(Conn) +# cellids, ncells = vertextocellmap(mesh) + +# # create the local shapes +# numverts = numvertices(mesh) + +# fns = Vector{Vector{Shape{T}}}() +# pos = Vector{P}() + +# sizehint!(fns, length(vertexlist)) +# sizehint!(pos, length(vertexlist)) +# for v in vertexlist + +# numshapes = ncells[v] +# numshapes == 0 && continue # skip detached vertices + +# shapes = Vector{Shape{T}}(undef,numshapes) +# for s in 1: numshapes +# c = cellids[v,s] +# cell = mesh.faces[c] +# if cell[1] == v +# shapes[s] = Shape(c, 1, T(1.0)) +# elseif cell[2] == v +# shapes[s] = Shape(c, 2, T(1.0)) +# else +# error("Junctions not supported") +# end +# end + +# push!(fns, shapes) +# push!(pos, mesh.vertices[v]) +# end - T = coordtype(mesh) - P = vertextype(mesh) - S = Shape{T} +# NF = 2 +# LagrangeBasis{1,0,NF}(geometry, fns, pos) +# end + +# function lagrangec0(mesh, vertexlist, ::Type{Val{2}}; order=1) + +# if order > 10 +# error("Maximal 1D Lagrange polynomial order supported is 10") +# end - fns = Vector{Vector{S}}() - pos = Vector{P}() - for (i,node) in enumerate(nodes) - fn = Vector{S}() - for k in nzrange(Conn,i) - cellid = rows[k] - refid = vals[k] - push!(fn, Shape(cellid, refid, T(1.0))) - end - push!(fns,fn) - push!(pos,cartesian(center(chart(nodes,node)))) - end +# T = coordtype(mesh) +# U = universedimension(mesh) +# P = vertextype(mesh) +# S = Shape{T} - Conn = connectivity(edges, mesh, abs) - rows = rowvals(Conn) - vals = nonzeros(Conn) +# geometry = mesh - for (i,edge) in enumerate(edges) - fn = Vector{S}() - for k in nzrange(Conn,i) - cellid = rows[k] - refid = vals[k] - push!(fn, Shape(cellid, 3+refid, T(1.0))) - end - push!(fns,fn) - push!(pos,cartesian(center(chart(edges,edge)))) - end +# cellids, ncells = vertextocellmap(mesh) + +# # create the local shapes +# numverts = numvertices(mesh) - NF = 6 - LagrangeBasis{2,0,NF}(mesh, fns, pos) -end +# fns = Vector{Vector{S}}() +# pos = Vector{P}() + +# sizehint!(fns, length(vertexlist)) +# sizehint!(pos, length(vertexlist)) +# for v in vertexlist + +# numshapes = ncells[v] +# numshapes == 0 && continue # skip detached vertices + +# shapes = Vector{S}(undef,numshapes) +# for s in 1: numshapes +# c = cellids[v,s] +# cell = mesh.faces[c] +# if cell[1] == v +# shapes[s] = Shape(c, 1, T(1.0)) +# elseif cell[2] == v +# shapes[s] = Shape(c, 2, T(1.0)) +# else +# error("Junctions not supported") +# end +# end + +# push!(fns, shapes) +# push!(pos, mesh.vertices[v]) +# end + +# NF = order + 1 + +# for (c,cell) in enumerate(mesh) +# ch = chart(mesh,cell) +# for r in 3:NF +# push!(fns, S[S(c,r,T(1.0))]) +# push!(pos, cartesian(neighborhood(ch, (r-2)/(NF-1)))) +# end +# end + +# LagrangeBasis{order,0,NF}(geometry, fns, pos) +# end + +# function lagrangec0d1(mesh, nodes::CompScienceMeshes.AbstractMesh{<:Any,1}) + +# Conn = connectivity(nodes, mesh, abs) +# rows = rowvals(Conn) +# vals = nonzeros(Conn) + +# T = coordtype(mesh) +# P = vertextype(mesh) +# S = Shape{T} + +# fns = Vector{Vector{S}}() +# pos = Vector{P}() +# for (i,node) in enumerate(nodes) +# fn = Vector{S}() +# for k in nzrange(Conn,i) +# cellid = rows[k] +# refid = vals[k] +# push!(fn, Shape(cellid, refid, T(1.0))) +# end +# push!(fns,fn) +# push!(pos,cartesian(center(chart(nodes,node)))) +# end + +# NF = dimension(mesh) + 1 + +# LagrangeBasis{1,0,NF}(mesh, fns, pos) +# end + +# function lagrangec0(mesh, nodes::CompScienceMeshes.AbstractMesh{<:Any,1}; order=1) + +# if order > 10 +# error("Maximal 1D Lagrange polynomial order supported is 10") +# end + +# Conn = connectivity(nodes, mesh, abs) +# rows = rowvals(Conn) +# vals = nonzeros(Conn) + +# T = coordtype(mesh) +# P = vertextype(mesh) +# S = Shape{T} + +# fns = Vector{Vector{S}}() +# pos = Vector{P}() + +# u = T(1.0) + +# for (i,node) in enumerate(nodes) +# fn = Vector{S}() +# for k in nzrange(Conn,i) +# cellid = rows[k] +# refid = vals[k] +# push!(fn, Shape(cellid, refid, u)) +# end +# push!(fns,fn) +# push!(pos,cartesian(center(chart(nodes,node)))) +# end + +# NF = order + 1 + +# for (c, cell) in enumerate(mesh) +# ch = chart(mesh, cell) +# for r in 3:NF +# push!(fns, S[S(c, r, u)]) +# push!(pos, cartesian(neighborhood(ch, SVector(u - (r-2)/(NF-1))))) +# end +# end + +# LagrangeBasis{order,0,NF}(mesh, fns, pos) +# end + +# function lagrangec0d2(mesh::CompScienceMeshes.AbstractMesh{U,3}, +# nodes::CompScienceMeshes.AbstractMesh{U,1}, +# edges::CompScienceMeshes.AbstractMesh{U,2}) where {U} + +# Conn = connectivity(nodes, mesh, abs) +# rows = rowvals(Conn) +# vals = nonzeros(Conn) + +# T = coordtype(mesh) +# P = vertextype(mesh) +# S = Shape{T} + +# fns = Vector{Vector{S}}() +# pos = Vector{P}() +# for (i,node) in enumerate(nodes) +# fn = Vector{S}() +# for k in nzrange(Conn,i) +# cellid = rows[k] +# refid = vals[k] +# push!(fn, Shape(cellid, refid, T(1.0))) +# end +# push!(fns,fn) +# push!(pos,cartesian(center(chart(nodes,node)))) +# end + +# Conn = connectivity(edges, mesh, abs) +# rows = rowvals(Conn) +# vals = nonzeros(Conn) + +# for (i,edge) in enumerate(edges) +# fn = Vector{S}() +# for k in nzrange(Conn,i) +# cellid = rows[k] +# refid = vals[k] +# push!(fn, Shape(cellid, 3+refid, T(1.0))) +# end +# push!(fns,fn) +# push!(pos,cartesian(center(chart(edges,edge)))) +# end + +# NF = 6 +# LagrangeBasis{2,0,NF}(mesh, fns, pos) +# end @@ -803,222 +1349,208 @@ function duallagrangec0d1(mesh, mesh2, pred, ::Type{Val{2}}) push!(pos, cartesian(center(chart(mesh, segment_coarse)))) end - NF = 2 - LagrangeBasis{1,0,NF}(geometry, fns, pos) -end - -""" - lagrangec0d2(mesh) -> basis - -Build lagrangec0d2 elements, including boundary vertices (meaning assuming a closed mesh). -""" - -function lagrangec0d2(mesh) - - #mark non connectet vertices - verts = skeleton(mesh, 0) - detached = trues(numvertices(mesh)) - for v in cells(verts) - detached[v.indices[1]] = false - end - - #identify - bnd = boundary(mesh) - bndverts = skeleton(bnd, 0) - notonbnd = trues(numvertices(mesh)) - for v in cells(bndverts) - notonbnd[v.indices[1]] = false - end - - vertexlist = findall(.!detached) - - edges = skeleton(mesh, 1; sort) - cps = cellpairs(mesh, edges, dropjunctionpair=true) - #only interior edges - ids = findall(x -> x>0, cps[2,:]) - - lagrangec0d2(mesh,vertexlist,cps,Val{dimension(mesh)+1}) -end - - -function lagrangec0d2(mesh, vertexlist::Vector, cellpairs::Array{Int,2}, ::Type{Val{3}}) - - T = coordtype(mesh) - U = universedimension(mesh) - - numvertices = length(vertexlist) - numedges = size(cellpairs,2) + NF = 2 + LagrangeBasis{1,0,NF}(geometry, fns, pos) +end - cellids, ncells = vertextocellmap(mesh) +# """ +# lagrangec0d2(mesh) -> basis - Cells = cells(mesh) - Verts = vertices(mesh) +# Build lagrangec0d2 elements, including boundary vertices (meaning assuming a closed mesh). +# """ - # create the local shapes - functions = Vector{Vector{Shape{Float64}}}(undef,numvertices+numedges) - positions = Vector{vertextype(mesh)}(undef,numvertices+numedges) - - #temp containers for vertex dof - fns = Vector{Vector{Shape{Float64}}}() - pos = Vector{vertextype(mesh)}() +# function lagrangec0d2(mesh) - sizehint!(fns, length(vertexlist)) - sizehint!(pos, length(vertexlist)) +# #mark non connectet vertices +# verts = skeleton(mesh, 0) +# detached = trues(numvertices(mesh)) +# for v in cells(verts) +# detached[v.indices[1]] = false +# end - #shape function associated with vetex dof - for v in vertexlist +# #identify +# bnd = boundary(mesh) +# bndverts = skeleton(bnd, 0) +# notonbnd = trues(numvertices(mesh)) +# for v in cells(bndverts) +# notonbnd[v.indices[1]] = false +# end - numshapes = ncells[v] - numshapes == 0 && continue +# vertexlist = findall(.!detached) + +# edges = skeleton(mesh, 1; sort) +# cps = cellpairs(mesh, edges, dropjunctionpair=true) +# #only interior edges +# ids = findall(x -> x>0, cps[2,:]) + +# lagrangec0d2(mesh,vertexlist,cps,Val{dimension(mesh)+1}) +# end - shapes = Vector{Shape{Float64}}(undef,numshapes) - for s in 1: numshapes - c = cellids[v,s] - # cell = mesh.faces[c] - cell = Cells[c] - localid = something(findfirst(isequal(v), cell.indices),0) - @assert localid != 0 +# function lagrangec0d2(mesh, vertexlist::Vector, cellpairs::Array{Int,2}, ::Type{Val{3}}) - shapes[s] = Shape(c, localid, 1.0) - end - #functions[v] = shapes - #positions[v] = Verts[v] +# T = coordtype(mesh) +# U = universedimension(mesh) - push!(fns, shapes) - push!(pos, mesh.vertices[v]) - end +# numvertices = length(vertexlist) +# numedges = size(cellpairs,2) - @assert length(fns) == numvertices +# cellids, ncells = vertextocellmap(mesh) - for i in 1:numvertices - functions[i] = fns[i] - positions[i] = pos[i] - end +# Cells = cells(mesh) +# Verts = vertices(mesh) - #shape function associated with edge dof - for i in 1:numedges - #internal edge - if cellpairs[2,i] > 0 - c1 = cellpairs[1,i]; cell1 = Cells[c1] #mesh.faces[c1] - c2 = cellpairs[2,i]; cell2 = Cells[c2] #mesh.faces[c2] - e1, e2 = getcommonedge(cell1.indices, cell2.indices) - functions[numvertices+i] = [ - Shape(c1, abs(e1)+3, 1.0), - Shape(c2, abs(e2)+3, 1.0)] - isct = intersect(cell1.indices, cell2.indices) - @assert length(isct) == 2 - @assert !(cell1.indices[abs(e1)] in isct) - @assert !(cell2.indices[abs(e2)] in isct) - - v1 = cell1.indices[mod1(abs(e1)+1,3)] - v2 = cell1.indices[mod1(abs(e1)+2,3)] - - edge = simplex(mesh.vertices[[v1,v2]]...) +# # create the local shapes +# functions = Vector{Vector{Shape{Float64}}}(undef,numvertices+numedges) +# positions = Vector{vertextype(mesh)}(undef,numvertices+numedges) + +# #temp containers for vertex dof +# fns = Vector{Vector{Shape{Float64}}}() +# pos = Vector{vertextype(mesh)}() + +# sizehint!(fns, length(vertexlist)) +# sizehint!(pos, length(vertexlist)) + +# #shape function associated with vetex dof +# for v in vertexlist + +# numshapes = ncells[v] +# numshapes == 0 && continue + +# shapes = Vector{Shape{Float64}}(undef,numshapes) +# for s in 1: numshapes +# c = cellids[v,s] +# # cell = mesh.faces[c] +# cell = Cells[c] + +# localid = something(findfirst(isequal(v), cell.indices),0) +# @assert localid != 0 + +# shapes[s] = Shape(c, localid, 1.0) +# end +# #functions[v] = shapes +# #positions[v] = Verts[v] + +# push!(fns, shapes) +# push!(pos, mesh.vertices[v]) +# end + +# @assert length(fns) == numvertices + +# for i in 1:numvertices +# functions[i] = fns[i] +# positions[i] = pos[i] +# end + +# #shape function associated with edge dof +# for i in 1:numedges +# #internal edge +# if cellpairs[2,i] > 0 +# c1 = cellpairs[1,i]; cell1 = Cells[c1] #mesh.faces[c1] +# c2 = cellpairs[2,i]; cell2 = Cells[c2] #mesh.faces[c2] +# e1, e2 = getcommonedge(cell1.indices, cell2.indices) +# functions[numvertices+i] = [ +# Shape(c1, abs(e1)+3, 1.0), +# Shape(c2, abs(e2)+3, 1.0)] +# isct = intersect(cell1.indices, cell2.indices) +# @assert length(isct) == 2 +# @assert !(cell1.indices[abs(e1)] in isct) +# @assert !(cell2.indices[abs(e2)] in isct) + +# v1 = cell1.indices[mod1(abs(e1)+1,3)] +# v2 = cell1.indices[mod1(abs(e1)+2,3)] + +# edge = simplex(mesh.vertices[[v1,v2]]...) - positions[numvertices+i] = cartesian(center(edge)) - #boundary edge - else - c1 = cellpairs[1,i]; cell1 = Cells[c1] - e1 = cellpairs[2,i] - functions[numvertices+i] = [ - Shape(c1, abs(e1)+3, +1.0)] - - v1 = cell1.indices[mod1(abs(e1)+1,3)] - v2 = cell1.indices[mod1(abs(e1)+2,3)] - edge = simplex(mesh.vertices[[v1,v2]]...) +# positions[numvertices+i] = cartesian(center(edge)) +# #boundary edge +# else +# c1 = cellpairs[1,i]; cell1 = Cells[c1] +# e1 = cellpairs[2,i] +# functions[numvertices+i] = [ +# Shape(c1, abs(e1)+3, +1.0)] + +# v1 = cell1.indices[mod1(abs(e1)+1,3)] +# v2 = cell1.indices[mod1(abs(e1)+2,3)] +# edge = simplex(mesh.vertices[[v1,v2]]...) - positions[numvertices+i] = cartesian(center(edge)) - end - end +# positions[numvertices+i] = cartesian(center(edge)) +# end +# end - NF = 6 - LagrangeBasis{2,0,NF}(mesh, functions, positions) -end +# NF = 6 +# LagrangeBasis{2,0,NF}(mesh, functions, positions) +# end -gradient(space::LagrangeBasis{1,0}, geo, fns) = NDLCCBasis(geo, fns) +# gradient(space::LagrangeBasis{1,0}, geo, fns) = NDLCCBasis(geo, fns) # gradient(space::LagrangeBasis{1,0}, geo::CompScienceMeshes.AbstractMesh{U,3} where {U}, fns) = NDBasis(geo, fns) #curl(space::LagrangeBasis{1,0}, geo, fns) = RTBasis(geo, fns) #curl(space::LagrangeBasis{2,0}, geo, fns) = BDMBasis(geo, fns) -gradient(space::LagrangeBasis{1,0,<:CompScienceMeshes.AbstractMesh{<:Any,2}}, geo, fns) = - LagrangeBasis{0,-1,1}(geo, fns, space.pos) +# gradient(space::LagrangeBasis{1,0,<:CompScienceMeshes.AbstractMesh{<:Any,2}}, geo, fns) = +# LagrangeBasis{0,-1,1}(geo, fns, space.pos) -gradient(space::LagrangeBasis{1,0,<:CompScienceMeshes.AbstractMesh{<:Any,3}}, geo, fns) = - NDBasis(geo, fns, space.pos) +# gradient(space::LagrangeBasis{1,0,<:CompScienceMeshes.AbstractMesh{<:Any,3}}, geo, fns) = +# NDBasis(geo, fns, space.pos) #curl(space::LagrangeBasis{2,0},geo, fns) = BDMBasis(geo, fns) #curl(space::LagrangeBasis{1,0}, geo, fns) = RTBasis(geo, fns) -function curl(X::LagrangeBasis{D,C,M,T,NF,P} , geo, fns) where {D,C,M,T,NF,P} - GWPDivSpace{T,M,Vector{P}}(geo, fns, deepcopy(positions(X)),D-1) -end +# function curl(X::LagrangeBasis{D,C,M,T,NF,P} , geo, fns) where {D,C,M,T,NF,P} +# GWPDivSpace{T,M,Vector{P}}(geo, fns, deepcopy(positions(X)),D-1) +# end -@testitem "curl - global" begin - using LinearAlgebra - using CompScienceMeshes - const CSM = CompScienceMeshes +# @testitem "curl - global" begin +# using LinearAlgebra +# using CompScienceMeshes +# const CSM = CompScienceMeshes - T = Float64 +# T = Float64 - m = CSM.meshrectangle(1.0, 1.0, 0.5, 3) - X = BEAST.lagrangec0(m; order=1) - curlX = BEAST.curl(X) +# m = CSM.meshrectangle(1.0, 1.0, 0.5, 3) +# X = BEAST.lagrangec0(m; order=1) +# curlX = BEAST.curl(X) - x = BEAST.refspace(X) - curlx = BEAST.refspace(curlX) +# x = BEAST.refspace(X) +# curlx = BEAST.refspace(curlX) - 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) +# 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 = zeros(T,3) - ϕp = x(p) - for sh in X.fns[i] - sh.cellid == cellid || continue - r1 += sh.coeff * ϕp[sh.refid].curl - end +# u = (0.2341, 0.4312) +# p = neighborhood(ch, u) - r2 = zeros(T,3) - ϕp = curlx(p) - for sh in curlX.fns[i] - sh.cellid == cellid || continue - r2 += sh.coeff * ϕp[sh.refid].value - end - global err = max(err, norm(r1-r2)) - end - end +# r1 = zeros(T,3) +# ϕp = x(p) +# for sh in X.fns[i] +# sh.cellid == cellid || continue +# r1 += sh.coeff * ϕp[sh.refid].curl +# end - @test err < sqrt(eps(T)) -end +# r2 = zeros(T,3) +# ϕp = curlx(p) +# for sh in curlX.fns[i] +# sh.cellid == cellid || continue +# r2 += sh.coeff * ϕp[sh.refid].value +# end +# global err = max(err, norm(r1-r2)) +# end +# end -# -# Sclar trace for Laggrange element based spaces -# -function strace(X::LagrangeBasis{1,0}, geo, fns::Vector) - # dimension of the boundary - n = dimension(geo) - # degree of the space - d = 1 - # number of local shape functions - NF = binomial(n+d,n) +# @test err < sqrt(eps(T)) +# end - trpos = deepcopy(positions(X)) - LagrangeBasis{1,0,NF}(geo, fns, trpos) -end @@ -1165,266 +1697,260 @@ function dual0forms_body(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}, refd, bn LagrangeBasis{1,0,3}(refd, bfs, pos) end -function lagrangecx(mesh::CompScienceMeshes.AbstractMesh{<:Any,2}; order) - - T = coordtype(mesh) - - P = vertextype(mesh) - S = Shape{T} - - fns = Vector{Vector{S}}() - pos = Vector{P}() - - u = one(T) - for (c,cell) in enumerate(mesh) - ch = chart(mesh,cell) - - push!(fns, S[S(c,1,u)]) - push!(pos, cartesian(neighborhood(ch, T(1.0)))) - - push!(fns, S[S(c,2,u)]) - push!(pos, cartesian(neighborhood(ch, T(0.0)))) - end - - NF = order + 1 - - for (c,cell) in enumerate(mesh) - ch = chart(mesh,cell) - for r in 3:NF - push!(fns, S[S(c,r,u)]) - push!(pos, cartesian(neighborhood(ch, u - (r-2)/(NF-1)))) - end - end - - return LagrangeBasis{order,-1,NF}(mesh, fns, pos) -end - -function lagrangecx(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}; order) - - T = coordtype(mesh) - NF = binomial(2+order, 2) - P = vertextype(mesh) - S = Shape{T} - - fns = Vector{Vector{S}}(undef, length(mesh) * NF) - pos = Vector{P}(undef, length(mesh) * NF) - - idx = 1 - u = one(T) - for (c,cell) in enumerate(mesh) - ch = chart(mesh,cell) - for r in 1:NF - fns[idx] = S[S(c,r,u)] - pos[idx] = cartesian(center(ch)) - idx += 1 - end - end - - return LagrangeBasis{order,-1,NF}(mesh, fns, pos) -end - - - -""" - localindices(dof, chart, i) - -Returns a vector of indices into the vector of local shape functions that correspond to -global degrees of freedom supported on sub-entity `i`, where the type of entity (nodes, -edge, face) is encoded in the type of 'dof'. -""" -function localindices(dof::_LagrangeGlobalNodesDoFs, chart::CompScienceMeshes.Simplex, - localspace, i) - - d = dof.order - [div((d+1)*(d+2),2), d+1, 1][[i]] -end - -function localindices(dof::_LagrangeGlobalEdgeDoFs, chart::CompScienceMeshes.Simplex, - localspace, i) - - d = dof.order - lids, lid = Int[], 0 - if i == 1 - for i in 0:d for j in 0:d - k = d - i - j - k < 0 && continue - lid += 1 - i != 0 && continue - j in (0,d) && continue - push!(lids, lid) - end end - elseif i == 2 - for i in 0:d for j in 0:d - k = d - i - j - k < 0 && continue - lid += 1 - i in (0,d) && continue - j != 0 && continue - push!(lids, lid) - end end - elseif i ==3 - for i in 0:d for j in 0:d - k = d - i - j - k < 0 && continue - lid += 1 - j in (0,d) && continue - k != 0 && continue - push!(lids, lid) - end end - else - error("wrong local edge index") - end - return lids -end - -function localindices(dof::_LagrangeGlobalFaceDoFs, chart::CompScienceMeshes.Simplex, - localspace, i) - - d = dof.order - lids, lid = Int[], 0 - for i in 0:d, j in 0:d - k = d - i - j - k < 0 && continue - lid += 1 - (i == 0 || j == 0 || k == 0) && continue - # @show (i,j,k) - push!(lids, lid) - end - return lids -end - -function localindices(dof::_LagrangeGlobalNodesDoFs, chart::CompScienceMeshes.Simplex, - localspace::LagrangeRefSpace{<:Real,2}, i) - return [1,3,6][[i]] -end - -function localindices(dof::_LagrangeGlobalEdgeDoFs, chart::CompScienceMeshes.Simplex, - localspace::LagrangeRefSpace{<:Real,2}, i) - return [5,4,2][[i]] -end - -function localindices(dof::_LagrangeGlobalFaceDoFs, chart::CompScienceMeshes.Simplex, - localspace::LagrangeRefSpace{<:Real,2}, i) - return [] -end - -function localindices(dof::_LagrangeGlobalNodesDoFs, chart::CompScienceMeshes.Simplex, - localspace::LagrangeRefSpace{<:Real,1}, i) - return [i] -end - -function localindices(dof::_LagrangeGlobalEdgeDoFs, chart::CompScienceMeshes.Simplex, - localspace::LagrangeRefSpace{<:Real,1}, i) - return [] -end - -function localindices(dof::_LagrangeGlobalFaceDoFs, chart::CompScienceMeshes.Simplex, - localspace::LagrangeRefSpace{<:Real,1}, i) - return [] -end - - - -function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}; order) - - T = coordtype(mesh) - atol = sqrt(eps(T)) - - P = vertextype(mesh) - S = Shape{T} - F = Vector{S} - - verts = skeleton(mesh, 0) - edges = skeleton(mesh, 1) - - C02 = connectivity(mesh, verts, abs); R02 = rowvals(C02); V02 = nonzeros(C02) - C12 = connectivity(mesh, edges, abs); R12 = rowvals(C12); V12 = nonzeros(C12) - - ne = order-1 - nf = div((order-2)*(order-1), 2) - - nV = length(verts) - nE = length(edges) * ne - nF = length(mesh) * nf - N = nV + nE + nF - - localspace = LagrangeRefSpace{T,order,3,binomial(2+order,2)}() - dom = domain(chart(mesh, first(mesh))) - localdim = numfunctions(localspace, dom) +# function lagrangecx(mesh::CompScienceMeshes.AbstractMesh{<:Any,2}; order) + +# T = coordtype(mesh) + +# P = vertextype(mesh) +# S = Shape{T} + +# fns = Vector{Vector{S}}() +# pos = Vector{P}() + +# u = one(T) +# for (c,cell) in enumerate(mesh) +# ch = chart(mesh,cell) + +# push!(fns, S[S(c,1,u)]) +# push!(pos, cartesian(neighborhood(ch, T(1.0)))) + +# push!(fns, S[S(c,2,u)]) +# push!(pos, cartesian(neighborhood(ch, T(0.0)))) +# end + +# NF = order + 1 + +# for (c,cell) in enumerate(mesh) +# ch = chart(mesh,cell) +# for r in 3:NF +# push!(fns, S[S(c,r,u)]) +# push!(pos, cartesian(neighborhood(ch, u - (r-2)/(NF-1)))) +# end +# end + +# return LagrangeBasis{order,-1,NF}(mesh, fns, pos) +# end + +# function lagrangecx(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}; order) + +# T = coordtype(mesh) +# NF = binomial(2+order, 2) +# P = vertextype(mesh) +# S = Shape{T} + +# fns = Vector{Vector{S}}(undef, length(mesh) * NF) +# pos = Vector{P}(undef, length(mesh) * NF) + +# idx = 1 +# u = one(T) +# for (c,cell) in enumerate(mesh) +# ch = chart(mesh,cell) +# for r in 1:NF +# fns[idx] = S[S(c,r,u)] +# pos[idx] = cartesian(center(ch)) +# idx += 1 +# end +# end + +# return LagrangeBasis{order,-1,NF}(mesh, fns, pos) +# end + + + +# """ +# localindices(dof, chart, i) + +# Returns a vector of indices into the vector of local shape functions that correspond to +# global degrees of freedom supported on sub-entity `i`, where the type of entity (nodes, +# edge, face) is encoded in the type of 'dof'. +# """ +# function localindices(dof::_LagrangeGlobalNodesDoFs, chart::CompScienceMeshes.Simplex, +# localspace, i) + +# d = dof.order +# [div((d+1)*(d+2),2), d+1, 1][[i]] +# end + +# function localindices(dof::_LagrangeGlobalEdgeDoFs, chart::CompScienceMeshes.Simplex, +# localspace, i) + +# d = dof.order +# lids, lid = Int[], 0 +# if i == 1 +# for i in 0:d for j in 0:d +# k = d - i - j +# k < 0 && continue +# lid += 1 +# i != 0 && continue +# j in (0,d) && continue +# push!(lids, lid) +# end end +# elseif i == 2 +# for i in 0:d for j in 0:d +# k = d - i - j +# k < 0 && continue +# lid += 1 +# i in (0,d) && continue +# j != 0 && continue +# push!(lids, lid) +# end end +# elseif i ==3 +# for i in 0:d for j in 0:d +# k = d - i - j +# k < 0 && continue +# lid += 1 +# j in (0,d) && continue +# k != 0 && continue +# push!(lids, lid) +# end end +# else +# error("wrong local edge index") +# end +# return lids +# end + +# function localindices(dof::_LagrangeGlobalFaceDoFs, chart::CompScienceMeshes.Simplex, +# localspace, i) + +# d = dof.order +# lids, lid = Int[], 0 +# for i in 0:d, j in 0:d +# k = d - i - j +# k < 0 && continue +# lid += 1 +# (i == 0 || j == 0 || k == 0) && continue +# # @show (i,j,k) +# push!(lids, lid) +# end +# return lids +# end + +# function localindices(dof::_LagrangeGlobalNodesDoFs, chart::CompScienceMeshes.Simplex, +# localspace::LagrangeRefSpace{<:Real,2}, i) +# return [1,3,6][[i]] +# end + +# function localindices(dof::_LagrangeGlobalEdgeDoFs, chart::CompScienceMeshes.Simplex, +# localspace::LagrangeRefSpace{<:Real,2}, i) +# return [5,4,2][[i]] +# end + +# function localindices(dof::_LagrangeGlobalFaceDoFs, chart::CompScienceMeshes.Simplex, +# localspace::LagrangeRefSpace{<:Real,2}, i) +# return [] +# end + +# function localindices(dof::_LagrangeGlobalNodesDoFs, chart::CompScienceMeshes.Simplex, +# localspace::LagrangeRefSpace{<:Real,1}, i) +# return [i] +# end + +# function localindices(dof::_LagrangeGlobalEdgeDoFs, chart::CompScienceMeshes.Simplex, +# localspace::LagrangeRefSpace{<:Real,1}, i) +# return [] +# end + +# function localindices(dof::_LagrangeGlobalFaceDoFs, chart::CompScienceMeshes.Simplex, +# localspace::LagrangeRefSpace{<:Real,1}, i) +# return [] +# end + + + +# function lagrangec0(mesh::CompScienceMeshes.AbstractMesh{<:Any,3}; order) + +# T = coordtype(mesh) +# atol = sqrt(eps(T)) + +# P = vertextype(mesh) +# S = Shape{T} +# F = Vector{S} + +# verts = skeleton(mesh, 0) +# edges = skeleton(mesh, 1) + +# C02 = connectivity(mesh, verts, abs); R02 = rowvals(C02); V02 = nonzeros(C02) +# C12 = connectivity(mesh, edges, abs); R12 = rowvals(C12); V12 = nonzeros(C12) + +# ne = order-1 +# nf = div((order-2)*(order-1), 2) + +# nV = length(verts) +# nE = length(edges) * ne +# nF = length(mesh) * nf +# N = nV + nE + nF + +# localspace = LagrangeRefSpace{T,order,3,binomial(2+order,2)}() +# dom = domain(chart(mesh, first(mesh))) +# localdim = numfunctions(localspace, dom) - d = order - fns = [S[] for n in 1:(nV+nE+nF)] - pos = fill(point(0,0,0), nV+nE+nF) - for cell in mesh - cell_ch = chart(mesh, cell) - V = R02[nzrange(C02,cell)] - I = V02[nzrange(C02,cell)] - for (i,v) in zip(I, V) - vertex_ch = chart(verts, v) - gids = [v] - lids = localindices(_LagrangeGlobalNodesDoFs(d), cell_ch, localspace, i) - v = globaldofs(vertex_ch, cell_ch, localspace, _LagrangeGlobalNodesDoFs(d)) - α = v[lids, :] - β = inv(α') - for i in axes(β,1) - for j in axes(β,2) - isapprox(β[i,j], 0; atol) && continue - push!(fns[gids[i]], S(cell, lids[j], β[i,j])) - end - end - end - - order < 2 && continue - E = R12[nzrange(C12,cell)] - I = V12[nzrange(C12,cell)] - for (i,e) in zip(I, E) - edge_ch = chart(edges, e) - gids = nV + (e-1)*ne + 1: nV + e*ne - lids = localindices(_LagrangeGlobalEdgeDoFs(d), cell_ch, localspace, i) - @assert length(lids) == length(gids) - v = globaldofs(edge_ch, cell_ch, localspace, _LagrangeGlobalEdgeDoFs(d)) - α = v[lids, :] - β = inv(α') - for i in axes(β,1) - for j in axes(β,2) - isapprox(β[i,j], 0; atol) && continue - push!(fns[gids[i]], S(cell, lids[j], β[i,j])) - end - end - end - - order < 3 && continue - F = [cell] - for (q,f) in enumerate(F) - face_ch = chart(mesh, f) - gids = nV + nE + (f-1)*nf + 1 : nV + nE + f*nf - lids = localindices(_LagrangeGlobalFaceDoFs(d), cell_ch, localspace, 1) - v = globaldofs(face_ch, cell_ch, localspace, _LagrangeGlobalFaceDoFs(d)) - α = v[lids, :] - # @show α - β = inv(α') - for i in axes(β,1) - for j in axes(β,2) - isapprox(β[i,j], 0; atol) && continue - push!(fns[gids[i]], S(cell, lids[j], β[i,j])) - end - end - end - end - - for v in verts pos[v] = cartesian(center(chart(verts, v))) end - for e in edges pos[nV + (e-1)*ne + 1: nV + e*ne] .= Ref(cartesian(center(chart(edges, e)))) end - for f in mesh pos[nV + nE + (f-1)*nf + 1: nV + nE + f*nf] .= Ref(cartesian(center(chart(mesh, f)))) end - - return LagrangeBasis{order,0,localdim}(mesh, fns, pos) -end +# d = order +# fns = [S[] for n in 1:(nV+nE+nF)] +# pos = fill(point(0,0,0), nV+nE+nF) +# for cell in mesh +# cell_ch = chart(mesh, cell) +# V = R02[nzrange(C02,cell)] +# I = V02[nzrange(C02,cell)] +# for (i,v) in zip(I, V) +# vertex_ch = chart(verts, v) +# gids = [v] +# lids = localindices(_LagrangeGlobalNodesDoFs(d), cell_ch, localspace, i) +# v = globaldofs(vertex_ch, cell_ch, localspace, _LagrangeGlobalNodesDoFs(d)) +# α = v[lids, :] +# β = inv(α') +# for i in axes(β,1) +# for j in axes(β,2) +# isapprox(β[i,j], 0; atol) && continue +# push!(fns[gids[i]], S(cell, lids[j], β[i,j])) +# end +# end +# end + +# order < 2 && continue +# E = R12[nzrange(C12,cell)] +# I = V12[nzrange(C12,cell)] +# for (i,e) in zip(I, E) +# edge_ch = chart(edges, e) +# gids = nV + (e-1)*ne + 1: nV + e*ne +# lids = localindices(_LagrangeGlobalEdgeDoFs(d), cell_ch, localspace, i) +# @assert length(lids) == length(gids) +# v = globaldofs(edge_ch, cell_ch, localspace, _LagrangeGlobalEdgeDoFs(d)) +# α = v[lids, :] +# β = inv(α') +# for i in axes(β,1) +# for j in axes(β,2) +# isapprox(β[i,j], 0; atol) && continue +# push!(fns[gids[i]], S(cell, lids[j], β[i,j])) +# end +# end +# end + +# order < 3 && continue +# F = [cell] +# for (q,f) in enumerate(F) +# face_ch = chart(mesh, f) +# gids = nV + nE + (f-1)*nf + 1 : nV + nE + f*nf +# lids = localindices(_LagrangeGlobalFaceDoFs(d), cell_ch, localspace, 1) +# v = globaldofs(face_ch, cell_ch, localspace, _LagrangeGlobalFaceDoFs(d)) +# α = v[lids, :] +# # @show α +# β = inv(α') +# for i in axes(β,1) +# for j in axes(β,2) +# isapprox(β[i,j], 0; atol) && continue +# push!(fns[gids[i]], S(cell, lids[j], β[i,j])) +# end +# end +# end +# end + +# for v in verts pos[v] = cartesian(center(chart(verts, v))) end +# for e in edges pos[nV + (e-1)*ne + 1: nV + e*ne] .= Ref(cartesian(center(chart(edges, e)))) end +# for f in mesh pos[nV + nE + (f-1)*nf + 1: nV + nE + f*nf] .= Ref(cartesian(center(chart(mesh, f)))) end + +# return LagrangeBasis{order,0,localdim}(mesh, fns, pos) +# end -function _surface(X) - C = unitfunctioncxd0(X.geo) - G = assemble(Identity(),X,X) - u = Vector(assemble(Identity(),X,C)[1:end,1]) - dot(u,G\u) -end @testitem "duallagrangec0d1" begin using CompScienceMeshes Γ = meshcuboid(1.0,1.0,1.0,0.5) diff --git a/src/bases/local/gwpdivlocal.jl b/src/bases/local/gwpdivlocal.jl index 1c8d024b..5393a5dd 100644 --- a/src/bases/local/gwpdivlocal.jl +++ b/src/bases/local/gwpdivlocal.jl @@ -121,7 +121,7 @@ end D = 4 NF = binomial(2+D, 2) gwp = BEAST.GWPDivRefSpace{T,D}() - lgx = BEAST.LagrangeRefSpace{T,D,3,10}() + lgx = BEAST.LagrangeRefSpace{T,D,3,15}() ch = CSM.simplex( point(1,0,0), diff --git a/src/bases/local/laglocal.jl b/src/bases/local/laglocal.jl index 2b140c88..5f7e5aa3 100644 --- a/src/bases/local/laglocal.jl +++ b/src/bases/local/laglocal.jl @@ -14,9 +14,7 @@ numfunctions(s::LagrangeRefSpace{T,Dg}, ch::CompScienceMeshes.ReferenceSimplex{D valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} = T # Evaluate constant lagrange elements on anything -(ϕ::LagrangeRefSpace{T,0})(tp) where {T} = SVector(((value=one(T), derivative=zero(T)),)) -(ϕ::LagrangeRefSpace{T,0,3})(tp) where {T} = SVector(((value=one(T), derivative=zero(T)),)) - +(ϕ::LagrangeRefSpace{T,0,2})(tp) where {T} = SVector(((value=one(T), derivative=zero(T)),)) # Evaluate linear Lagrange elements on a segment # The derivative denotes the tangential derivative @@ -210,28 +208,189 @@ function (f::LagrangeRefSpace{T,10,2})(mp) where T ) end -# Evaluete linear lagrange elements on a triangle -function (f::LagrangeRefSpace{T,1,3})(t) where T - u,v,w, = barycentric(t) - j = jacobian(t) - p = t.patch - σ = sign(dot(normal(t), cross(p[1]-p[3],p[2]-p[3]))) - SVector( - (value=u, curl=σ*(p[3]-p[2])/j), - (value=v, curl=σ*(p[1]-p[3])/j), - (value=w, curl=σ*(p[2]-p[1])/j)) + + +# Evaluate Lagrange element on a triangle with the surface curl, for arbitrary degree +@generated function (::LagrangeRefSpace{T,Degree,3,NF})(mp) where {T, Degree, NF} + + u = :(parametric(mp)[1]) + v = :(parametric(mp)[2]) + w = :(1 - $u - $v) + + tu = :(tangents(mp,1)) + tv = :(tangents(mp,2)) + jd = :(jacobian(mp)) + + nodes = :() + for i in 0:Degree + nodes = :($nodes..., T($i/($Degree))) + end + + vals = :() + diffs = :() + for k in 1:Degree+1 + for j in 1:Degree+1 + for i in 1:Degree+1 + i + j + k == Degree+3 || continue + + pp = gen_lagpoly(nodes,i,u,1,i,T) + pq = gen_lagpoly(nodes,j,v,1,j,T) + pr = gen_lagpoly(nodes,k,w,1,k,T) + val = :($pp * $pq * $pr) + + vals = :($vals..., $val) + + diffu = :(zero(T)) + diffv = :(zero(T)) + + pl = gen_lagpoly_diff(nodes,i,u,1,i,T) + diffu = :($pl*$pq*$pr) + + pm = gen_lagpoly_diff(nodes,j,v,1,j,T) + diffv = :($pm*$pp*$pr) + + pn = gen_lagpoly_diff(nodes,k,w,1,k,T) + diffu = :($diffu - $pn*$pp*$pq) + diffv = :($diffv - $pn*$pp*$pq) + + diffs = :($diffs..., ($diffu, $diffv)) + end + end + end + + ex = :(SVector{NF}(())) + + for i in 1:NF + push!(ex.args[2].args, :(value=$vals[$i], curl=($diffs[$i][2]*$tu-$diffs[$i][1]*$tv)/$jd)) + end + + return ex +end + + +""" + localindices(dof, chart, i) + +Returns a vector of indices into the vector of local shape functions that correspond to +global degrees of freedom supported on sub-entity `i`, where the type of entity (nodes, +edge, face) is encoded in the type of 'dof'. +""" +function localindices(dof::_LagrangeGlobalNodesDoFs, chart::CompScienceMeshes.Simplex{<:Any, 2}, + localspace, i) + + d = dof.order + [1, d+1, div((d+1)*(d+2),2)][[i]] +end + +function localindices(dof::_LagrangeGlobalEdgeDoFs, chart::CompScienceMeshes.Simplex{<:Any, 2}, + localspace, i) + + d = dof.order + lids, lid = Int[], 0 + if i == 3 + for i in 0:d for j in 0:d + k = d - i - j + k < 0 && continue + lid += 1 + i != 0 && continue + j in (0,d) && continue + push!(lids, lid) + end end + elseif i == 2 + for i in 0:d for j in 0:d + k = d - i - j + k < 0 && continue + lid += 1 + i in (0,d) && continue + j != 0 && continue + push!(lids, lid) + end end + elseif i == 1 + for i in 0:d for j in 0:d + k = d - i - j + k < 0 && continue + lid += 1 + j in (0,d) && continue + k != 0 && continue + push!(lids, lid) + end end + else + error("wrong local edge index") + end + return lids end +function localindices(dof::_LagrangeGlobalFaceDoFs, chart::CompScienceMeshes.Simplex{<:Any, 2}, + localspace, i) + d = dof.order + lids, lid = Int[], 0 + for i in 0:d, j in 0:d + k = d - i - j + k < 0 && continue + lid += 1 + (i == 0 || j == 0 || k == 0) && continue + # @show (i,j,k) + push!(lids, lid) + end + return lids +end -# Evaluate constant Lagrange elements on a triangle, with their curls -function (f::LagrangeRefSpace{T,0,3})(t, ::Type{Val{:withcurl}}) where T - i = one(T) - z = zero(cartesian(t)) - SVector(((value=i, curl=z,),)) + +function (ϕ::LagrangeRefSpace{T, 1, 4, 4})(lag) where T + + u, v, w = parametric(lag) + + tu = tangents(lag, 1) + tv = tangents(lag, 2) + tw = tangents(lag, 3) + + B = [tu tv tw] + A = inv(transpose(B)) + + # gradient in u,v,w (unit tetrahedron) + gr1=SVector{3, T}(1.0, 0.0, 0.0) + gr2=SVector{3, T}(0.0, 1.0, 0.0) + gr3=SVector{3, T}(0.0, 0.0, 1.0) + gr4=SVector{3, T}(-1.0, -1.0, -1.0) + + return SVector(( + (value = u, gradient = A*gr1), + (value = v, gradient = A*gr2), + (value = w, gradient = A*gr3), + (value = T(1.0)-u-v-w, gradient = A*gr4) + )) +end + +function localindices(dof::_LagrangeGlobalNodesDoFs, chart::CompScienceMeshes.Simplex{<:Any, 3}, + localspace, i) + [1, 2, 3, 4][[i]] end + +# Evaluete linear lagrange elements on a triangle +# function (f::LagrangeRefSpace{T,1,3})(t) where T +# u,v,w, = barycentric(t) + +# j = jacobian(t) +# p = t.patch +# σ = sign(dot(normal(t), cross(p[1]-p[3],p[2]-p[3]))) +# SVector( +# (value=u, curl=σ*(p[3]-p[2])/j), +# (value=v, curl=σ*(p[1]-p[3])/j), +# (value=w, curl=σ*(p[2]-p[1])/j)) +# end + + + +# Evaluate constant Lagrange elements on a triangle, with their curls +# function (f::LagrangeRefSpace{T,0,3})(t, ::Type{Val{:withcurl}}) where T +# i = one(T) +# z = zero(cartesian(t)) +# SVector(((value=i, curl=z,),)) +# end + #= Replaced by generalized curl using GWP functions function curl(ref::LagrangeRefSpace{T,1,3} where {T}, sh, el) sh1 = Shape(sh.cellid, mod1(sh.refid+1,3), -sh.coeff) @@ -298,6 +457,69 @@ function gradient(ref::LagrangeRefSpace{T,1,2}, sh, seg) where {T} end +function curl(localspace::LagrangeRefSpace, sh, ch) + fns = curl(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 curl(localspace::LagrangeRefSpace{T,D,3,N}, cellid::Int, ch) where {T,D,N} + function fields(p) + map(localspace(p)) do x + x.curl + end + end + atol = sqrt(eps(T)) + gwp = GWPDivRefSpace{T,D-1}() + coeffs = interpolate(fields, gwp, 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(coeffs[r,c]) > atol) + end + return SVector{N}(A) +end + + + +@testitem "curl - chartwise" begin + using LinearAlgebra + using CompScienceMeshes + const CSM = CompScienceMeshes + + T = Float64 + D = 1 + NF = binomial(2+D, 2) + gwp = BEAST.GWPDivRefSpace{T,D-1}() + lgc = BEAST.LagrangeRefSpace{T,D,3,NF}() + + ch = CSM.simplex( + point(1,0,0), + point(0,1,0), + point(0,0,0)) + + curlfns = BEAST.curl(lgc, 1, ch) + + p = neighborhood(ch, (0.2341, 0.4321)) + gwp_vals = gwp(p) + lgc_vals = lgc(p) + + err = similar(Vector{T}, axes(gwp_vals)) + for i in eachindex(lgc_vals) + val1 = lgc_vals[i].curl + val2 = zero(val1) + for sh in curlfns[i] + val2 += sh.coeff * gwp_vals[sh.refid].value + end + err[i] = norm(val1 - val2) + end + atol = sqrt(eps(T)) + @test all(err .< atol) +end + + function strace(x::LagrangeRefSpace, cell, localid, face) @@ -372,109 +594,6 @@ function restrict(f::LagrangeRefSpace{T,1}, dom1, dom2) where T end - -## Quadratic Lagrange element on a triangle -function (f::LagrangeRefSpace{T,2,3})(t) where T - u,v,w, = barycentric(t) - - j = jacobian(t) - p = t.patch - - - σ = sign(dot(normal(t), cross(p[1]-p[3],p[2]-p[3]))) - SVector( - (value=u*(2*u-1), curl=σ*(p[3]-p[2])*(4u-1)/j), - (value=4*u*v, curl=4*σ*(u*(p[1]-p[3])+v*(p[3]-p[2]))/j), - (value=v*(2*v-1), curl=σ*(p[1]-p[3])*(4v-1)/j), - (value=4*w*u, curl=4*σ*(w*(p[3]-p[2])+u*(p[2]-p[1]))/j), - (value=4*v*w, curl=4*σ*(w*(p[1]-p[3])+v*(p[2]-p[1]))/j), - (value=w*(2*w-1), curl=σ*(p[2]-p[1])*(4w-1)/j), - ) -end - -#= Not sure if that is correct anymore -function curl(ref::LagrangeRefSpace{T,2,3} where {T}, sh, el) - #curl of lagc0d2 as combination of bdm functions - z=zero(typeof(sh.coeff)) - if sh.refid < 4 - sh1 = Shape(sh.cellid, mod1(2*sh.refid+1,6), +sh.coeff) - sh2 = Shape(sh.cellid, mod1(2*sh.refid+2,6), -3*sh.coeff) - sh3 = Shape(sh.cellid, mod1(2*sh.refid+3,6), +3*sh.coeff) - sh4 = Shape(sh.cellid, mod1(2*sh.refid+4,6), -sh.coeff) - else - sh1 = Shape(sh.cellid, mod1(2*sh.refid+4,6), z*sh.coeff) - sh2 = Shape(sh.cellid, mod1(2*sh.refid+5,6), -4*sh.coeff) - sh3 = Shape(sh.cellid, mod1(2*sh.refid+6,6), +4*sh.coeff) - sh4 = Shape(sh.cellid, mod1(2*sh.refid+7,6), z*sh.coeff) - end - return [sh1, sh2, sh3, sh4] -end -=# - - -function curl(localspace::LagrangeRefSpace, sh, ch) - fns = curl(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 curl(localspace::LagrangeRefSpace{T,D,3,N}, cellid::Int, ch) where {T,D,N} - function fields(p) - map(localspace(p)) do x - x.curl - end - end - atol = sqrt(eps(T)) - gwp = GWPDivRefSpace{T,D-1}() - coeffs = interpolate(fields, gwp, 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(coeffs[r,c]) > atol) - end - return SVector{N}(A) -end - - - -@testitem "curl - chartwise" begin - using LinearAlgebra - using CompScienceMeshes - const CSM = CompScienceMeshes - - T = Float64 - D = 1 - NF = binomial(2+D, 2) - gwp = BEAST.GWPDivRefSpace{T,D-1}() - lgc = BEAST.LagrangeRefSpace{T,D,3,NF}() - - ch = CSM.simplex( - point(1,0,0), - point(0,1,0), - point(0,0,0)) - - curlfns = BEAST.curl(lgc, 1, ch) - - p = neighborhood(ch, (0.2341, 0.4321)) - gwp_vals = gwp(p) - lgc_vals = lgc(p) - - err = similar(Vector{T}, axes(gwp_vals)) - for i in eachindex(lgc_vals) - val1 = lgc_vals[i].curl - val2 = zero(val1) - for sh in curlfns[i] - val2 += sh.coeff * gwp_vals[sh.refid].value - end - err[i] = norm(val1 - val2) - end - atol = sqrt(eps(T)) - @test all(err .< atol) -end - - function restrict(f::LagrangeRefSpace{T,2}, dom1, dom2) where T D = numfunctions(f) @@ -520,177 +639,197 @@ function restrict(f::LagrangeRefSpace{T,2}, dom1, dom2) where T end -const _vert_perms_lag = [ - (1,2,3), - (2,3,1), - (3,1,2), - (2,1,3), - (1,3,2), - (3,2,1), -] - -const _dof_perms_lag0 = [ - (1), - (1), - (1), - (1), - (1), - (1), -] -const _dof_perms_lag1 = [ - (1,2,3), - (3,1,2), - (2,3,1), - (2,1,3), - (1,3,2), - (3,2,1), -] - -function dof_permutation(::LagrangeRefSpace{<:Any,0}, vert_permutation) - i = findfirst(==(tuple(vert_permutation...)), _vert_perms_lag) - return _dof_perms_lag0[i] -end - -function dof_permutation(::LagrangeRefSpace{<:Any,1}, vert_permutation) - i = findfirst(==(tuple(vert_permutation...)), _vert_perms_lag) - return _dof_perms_lag1[i] -end - -function dof_perm_matrix(::LagrangeRefSpace{<:Any,0}, vert_permutation) - i = findfirst(==(tuple(vert_permutation...)), _vert_perms_rt) - @assert i != nothing - return _dof_lag0perm_matrix[i] -end - -function dof_perm_matrix(::LagrangeRefSpace{<:Any,1}, vert_permutation) - i = findfirst(==(tuple(vert_permutation...)), _vert_perms_rt) - @assert i != nothing - return _dof_rtperm_matrix[i] -end - -const _dof_lag0perm_matrix = [ - @SMatrix[1], # 1. {1,2,3} - @SMatrix[1], # 2. {2,3,1} - @SMatrix[1], # 3. {3,1,2} - @SMatrix[1], # 4. {2,1,3} - @SMatrix[1], # 5. {1,3,2} - @SMatrix[1] # 6. {3,2,1} -] -function (ϕ::LagrangeRefSpace{T, 1, 4, 4})(lag) where T +# Quadratic Lagrange element on a triangle +# function (f::LagrangeRefSpace{T,2,3})(t) where T +# u,v,w, = barycentric(t) - u, v, w = parametric(lag) +# j = jacobian(t) +# p = t.patch - tu = tangents(lag, 1) - tv = tangents(lag, 2) - tw = tangents(lag, 3) - - B = [tu tv tw] - A = inv(transpose(B)) - # gradient in u,v,w (unit tetrahedron) - gr1=SVector{3, T}(1.0, 0.0, 0.0) - gr2=SVector{3, T}(0.0, 1.0, 0.0) - gr3=SVector{3, T}(0.0, 0.0, 1.0) - gr4=SVector{3, T}(-1.0, -1.0, -1.0) +# σ = sign(dot(normal(t), cross(p[1]-p[3],p[2]-p[3]))) +# SVector( +# (value=u*(2*u-1), curl=σ*(p[3]-p[2])*(4u-1)/j), +# (value=4*u*v, curl=4*σ*(u*(p[1]-p[3])+v*(p[3]-p[2]))/j), +# (value=v*(2*v-1), curl=σ*(p[1]-p[3])*(4v-1)/j), +# (value=4*w*u, curl=4*σ*(w*(p[3]-p[2])+u*(p[2]-p[1]))/j), +# (value=4*v*w, curl=4*σ*(w*(p[1]-p[3])+v*(p[2]-p[1]))/j), +# (value=w*(2*w-1), curl=σ*(p[2]-p[1])*(4w-1)/j), +# ) +# end - return SVector(( - (value = u, gradient = A*gr1), - (value = v, gradient = A*gr2), - (value = w, gradient = A*gr3), - (value = T(1.0)-u-v-w, gradient = A*gr4) - )) +#= Not sure if that is correct anymore +function curl(ref::LagrangeRefSpace{T,2,3} where {T}, sh, el) + #curl of lagc0d2 as combination of bdm functions + z=zero(typeof(sh.coeff)) + if sh.refid < 4 + sh1 = Shape(sh.cellid, mod1(2*sh.refid+1,6), +sh.coeff) + sh2 = Shape(sh.cellid, mod1(2*sh.refid+2,6), -3*sh.coeff) + sh3 = Shape(sh.cellid, mod1(2*sh.refid+3,6), +3*sh.coeff) + sh4 = Shape(sh.cellid, mod1(2*sh.refid+4,6), -sh.coeff) + else + sh1 = Shape(sh.cellid, mod1(2*sh.refid+4,6), z*sh.coeff) + sh2 = Shape(sh.cellid, mod1(2*sh.refid+5,6), -4*sh.coeff) + sh3 = Shape(sh.cellid, mod1(2*sh.refid+6,6), +4*sh.coeff) + sh4 = Shape(sh.cellid, mod1(2*sh.refid+7,6), z*sh.coeff) + end + return [sh1, sh2, sh3, sh4] end +=# -# Evaluate higher order Lagrange elements on triangles -# TODO: Optimise using code generation -function (ϕ::LagrangeRefSpace{T,Degree,3})(p) where {T,Degree} - u, v = parametric(p) - w = 1 - u - v - idx = 0 - suppdim = 2 - localdim = binomial(suppdim+Degree, suppdim) - vals = T[] - diffus = T[] - diffvs = T[] - D1 = Degree + 1 - s = range(zero(T), one(T), length=D1) - for i in 0:Degree - ui = i/Degree - for j in 0:Degree - vj = j/Degree - for k in 0:Degree - wk = k/Degree - i + j + k == Degree || continue - - prod_p = one(T) - for p in 0:i-1 - up = p / Degree - prod_p *= (u-up) / (ui-up) - end - prod_q = one(T) - for q in 0:j-1 - vq = q / Degree - prod_q *= (v-vq) / (vj-vq) - end - prod_r = one(T) - for r in 0:k-1 - wr = r / Degree - prod_r *= (w-wr) / (wk-wr) - end - push!(vals, prod_p * prod_q * prod_r) - - diffu = zero(T) - diffv = zero(T) - for l in 0:i-1 - ul = l/Degree - prod_pl = one(T) - for p in 0:i-1 - p == l && continue - up = p/Degree - prod_pl *= (u-up) / (ui-up) - end - diffu += prod_pl * prod_q * prod_r / (ui-ul) - end - for m in 0:j-1 - vm = m/Degree - prod_qm = one(T) - for q in 0:j-1 - q == m && continue - vq = q/Degree - prod_qm *= (v-vq) / (vj-vq) - end - diffv += prod_p * prod_qm * prod_r / (vj-vm) - end - for n in 0:k-1 - wn = n/Degree - prod_rn = one(T) - for r in 0:k-1 - r == n && continue - wr = r/Degree - prod_rn *= (w-wr) / (wk-wr) - end - diffu -= prod_p * prod_q * prod_rn / (wk-wn) - diffv -= prod_p * prod_q * prod_rn / (wk-wn) - end - - push!(diffus, diffu) - push!(diffvs, diffv) - - idx += 1 - end end end +# const _vert_perms_lag = [ +# (1,2,3), +# (2,3,1), +# (3,1,2), +# (2,1,3), +# (1,3,2), +# (3,2,1), +# ] + +# const _dof_perms_lag0 = [ +# (1), +# (1), +# (1), +# (1), +# (1), +# (1), +# ] +# const _dof_perms_lag1 = [ +# (1,2,3), +# (3,1,2), +# (2,3,1), +# (2,1,3), +# (1,3,2), +# (3,2,1), +# ] + +# function dof_permutation(::LagrangeRefSpace{<:Any,0}, vert_permutation) +# i = findfirst(==(tuple(vert_permutation...)), _vert_perms_lag) +# return _dof_perms_lag0[i] +# end + +# function dof_permutation(::LagrangeRefSpace{<:Any,1}, vert_permutation) +# i = findfirst(==(tuple(vert_permutation...)), _vert_perms_lag) +# return _dof_perms_lag1[i] +# end + +# function dof_perm_matrix(::LagrangeRefSpace{<:Any,0}, vert_permutation) +# i = findfirst(==(tuple(vert_permutation...)), _vert_perms_rt) +# @assert i != nothing +# return _dof_lag0perm_matrix[i] +# end + +# function dof_perm_matrix(::LagrangeRefSpace{<:Any,1}, vert_permutation) +# i = findfirst(==(tuple(vert_permutation...)), _vert_perms_rt) +# @assert i != nothing +# return _dof_rtperm_matrix[i] +# end + +# const _dof_lag0perm_matrix = [ +# @SMatrix[1], # 1. {1,2,3} +# @SMatrix[1], # 2. {2,3,1} +# @SMatrix[1], # 3. {3,1,2} +# @SMatrix[1], # 4. {2,1,3} +# @SMatrix[1], # 5. {1,3,2} +# @SMatrix[1] # 6. {3,2,1} +# ] + + + + +# Evaluate higher order Lagrange elements on triangles +# # TODO: Optimise using code generation +# function (ϕ::LagrangeRefSpaceRef{T,Degree,3})(p) where {T,Degree} + +# u, v = parametric(p) +# w = 1 - u - v +# idx = 0 + +# suppdim = 2 +# localdim = binomial(suppdim+Degree, suppdim) +# vals = T[] +# diffus = T[] +# diffvs = T[] + +# D1 = Degree + 1 +# s = range(zero(T), one(T), length=D1) +# for i in 0:Degree +# ui = i/Degree +# for j in 0:Degree +# vj = j/Degree +# for k in 0:Degree +# wk = k/Degree +# i + j + k == Degree || continue + +# prod_p = one(T) +# for p in 0:i-1 +# up = p / Degree +# prod_p *= (u-up) / (ui-up) +# end +# prod_q = one(T) +# for q in 0:j-1 +# vq = q / Degree +# prod_q *= (v-vq) / (vj-vq) +# end +# prod_r = one(T) +# for r in 0:k-1 +# wr = r / Degree +# prod_r *= (w-wr) / (wk-wr) +# end +# push!(vals, prod_p * prod_q * prod_r) + +# diffu = zero(T) +# diffv = zero(T) +# for l in 0:i-1 +# ul = l/Degree +# prod_pl = one(T) +# for p in 0:i-1 +# p == l && continue +# up = p/Degree +# prod_pl *= (u-up) / (ui-up) +# end +# diffu += prod_pl * prod_q * prod_r / (ui-ul) +# end +# for m in 0:j-1 +# vm = m/Degree +# prod_qm = one(T) +# for q in 0:j-1 +# q == m && continue +# vq = q/Degree +# prod_qm *= (v-vq) / (vj-vq) +# end +# diffv += prod_p * prod_qm * prod_r / (vj-vm) +# end +# for n in 0:k-1 +# wn = n/Degree +# prod_rn = one(T) +# for r in 0:k-1 +# r == n && continue +# wr = r/Degree +# prod_rn *= (w-wr) / (wk-wr) +# end +# diffu -= prod_p * prod_q * prod_rn / (wk-wn) +# diffv -= prod_p * prod_q * prod_rn / (wk-wn) +# end + +# push!(diffus, diffu) +# push!(diffvs, diffv) + +# idx += 1 +# end end end - tu = tangents(p,1) - tv = tangents(p,2) - j = jacobian(p) - NF = length(vals) - SVector{NF}([(value=f, curl=(-dv*tu+du*tv)/j) for (f,du,dv) in zip(vals, diffus, diffvs)]) -end +# tu = tangents(p,1) +# tv = tangents(p,2) +# j = jacobian(p) +# NF = length(vals) +# SVector{NF}([(value=f, curl=(-dv*tu+du*tv)/j) for (f,du,dv) in zip(vals, diffus, diffvs)]) +# end # fields[i] ≈ sum(Q[j,i] * interpolant[j].value for j in 1:numfunctions(interpolant)) function interpolate(fields, interpolant::LagrangeRefSpace{T,Degree,3}, chart) where {T,Degree} @@ -701,9 +840,9 @@ function interpolate(fields, interpolant::LagrangeRefSpace{T,Degree,3}, chart) w I = 0:Degree s = range(0,1,length=Degree+1) Is = zip(I,s) - for (i,ui) in Is + for (k,wk) in Is for (j,vj) in Is - for (k,wk) in Is + for (i,ui) in Is i + j + k == Degree || continue @assert ui + vj + wk ≈ 1 p = neighborhood(chart, (ui,vj)) @@ -721,41 +860,42 @@ function interpolate(fields, interpolant::LagrangeRefSpace{T,Degree,3}, chart) w return Q end -function interpolate(fields, interpolant::LagrangeRefSpace{T,1,3}, chart) where {T} - vals = Vector{Vector{T}}() +# function interpolate(fields, interpolant::LagrangeRefSpace{T,1,3}, chart) where {T} +# vals = Vector{Vector{T}}() - p1 = neighborhood(chart, (1,0)) - p2 = neighborhood(chart, (0,1)) - p3 = neighborhood(chart, (0,0)) - push!(vals, fields(p1)) - push!(vals, fields(p2)) - push!(vals, fields(p3)) +# p1 = neighborhood(chart, (1,0)) +# p2 = neighborhood(chart, (0,1)) +# p3 = neighborhood(chart, (0,0)) +# push!(vals, fields(p1)) +# push!(vals, fields(p2)) +# push!(vals, fields(p3)) - # 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 +# # 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 function interpolate!(out, fields, interpolant::LagrangeRefSpace{T,Degree,3}, chart) where {T,Degree} Is = zip((0:Degree), range(0,1,length=Degree+1)) idx = 0 - for (i,ui) in Is + for (k,wk) in Is for (j,vj) in Is - for (k,wk) in Is + for (i,ui) in Is i + j + k == Degree || continue; idx += 1 @assert ui + vj + wk ≈ 1 p = neighborhood(chart, (ui,vj)) vals = fields(p) for (g, val) in zip(axes(out, 1), vals) out[g,idx] = val -end end end end end + end end end end +end function interpolate!(out, fields, interpolant::LagrangeRefSpace{T,0,3}, chart) where {T} p = center(chart) diff --git a/test/test_basis.jl b/test/test_basis.jl index 90358764..15be4738 100644 --- a/test/test_basis.jl +++ b/test/test_basis.jl @@ -10,7 +10,7 @@ for T in [Float32, Float64] κ = ω = T(1.0) Γ = meshsegment(T(1.0), T(0.5)) - X = lagrangec0d1(Γ) + X = lagrangec0d1(Γ;dirichlet=true) @test numvertices(Γ)-2 == numfunctions(X) hypersingular = Helmholtz2D.hypersingular(; wavenumber=κ) @@ -75,7 +75,7 @@ using Test for T in [Float32, Float64] m = meshrectangle(T(1.0), T(1.0), T(0.5), 3) - X = lagrangec0d1(m) + X = lagrangec0d1(m,dirichlet=true) x = refspace(X) dom = domain(chart(m, first(m))) @@ -127,7 +127,7 @@ for T in [Float64] m = Mesh(V,C) b = Mesh([p1,p2], [CompScienceMeshes.SimplexGraph(1,2)]) - X = lagrangec0d1(m, boundary(m)) + X = lagrangec0d1(m) @test numfunctions(X) == 3 Y = BEAST.strace(X, b) @@ -161,7 +161,7 @@ for T in [Float64] m3 = CompScienceMeshes.rotate(m1, T(1.0π)*[1,0,0]) m = weld(m1, m2, m3) - X = lagrangec0d1(m) + X = lagrangec0d1(m,dirichlet=true) x = refspace(X) @test numfunctions(X) == 1 @test length(X.fns[1]) == 9 @@ -233,10 +233,10 @@ end ## Test gradient and curl of continuous lagrange elements m = meshrectangle(1.0, 1.0, 0.5, 3) -int_nodes = submesh(!in(skeleton(boundary(m),0)), skeleton(m,0)) -@test length(int_nodes) == 1 +#int_nodes = submesh(!in(skeleton(boundary(m),0)), skeleton(m,0)) +#@test length(int_nodes) == 1 -lag = lagrangec0d1(m, int_nodes) +lag = lagrangec0d1(m, dirichlet=true) @test numfunctions(lag) == 1 rs = refspace(lag) @@ -269,7 +269,7 @@ m = Mesh([ point(0,0,0)], [CompScienceMeshes.SimplexGraph(1,2,3,4)]) -lag = lagrangec0d1(m, skeleton(m,0)) +lag = lagrangec0d1(m) @test numfunctions(lag) == 4 gradlag = gradient(lag) diff --git a/test/test_gridfunction.jl b/test/test_gridfunction.jl index d23e2ae8..669a26b7 100644 --- a/test/test_gridfunction.jl +++ b/test/test_gridfunction.jl @@ -11,13 +11,13 @@ a = U(1) # Cube between (3,3,3) and (4,4,4) Γ = translate( - CompScienceMeshes.meshcuboid(a,a,a,U(1.0); generator=:gmsh), + CompScienceMeshes.meshcuboid(a,a,a,U(0.5)), SVector(3.0,3.0,3.0)) C0 = lagrangec0d1(Γ) # We interpolate the function f(x, y, z) = z -coeffs = [v[3] for v in vertices(Γ)] +coeffs = [v[3] for v in C0.pos] gfc0 = BEAST.FEMFunction(coeffs, C0) @@ -50,7 +50,7 @@ constgfcx = BEAST.FEMFunction([2.0 for c in cells(Γ)], CX) # We interpolate the function f(x, y, z) = z # with piecewise constant lagrange elements -coeffsCx = [((Γ.vertices[c[1]] + Γ.vertices[c[2]] + Γ.vertices[c[3]])/3)[3] for c in cells(Γ)] +coeffsCx = [v[3] for v in CX.pos] #.+ [((Γ.vertices[c[1]] + Γ.vertices[c[2]] + Γ.vertices[c[3]])/3)[3] for c in cells(Γ)] gfcx = BEAST.FEMFunction(coeffsCx, CX) @@ -61,7 +61,7 @@ gfcx = BEAST.FEMFunction(coeffsCx, CX) #difflincomb1 = BEAST.LinearCombinationOfAbstractMeshFunctions([1.0; -1.0], [gfcx, gfc0]) difflincomb1 = gfcx - gfc0 -@test BEAST.Lp_integrate(difflincomb1)/BEAST.Lp_integrate(gfc0) ≈ 0.03866223365135301 +@test BEAST.Lp_integrate(difflincomb1)/BEAST.Lp_integrate(gfc0) ≈ 0.02733832759069041 ## Test GlobalFunction f(x) = x[3] @@ -83,7 +83,7 @@ idxfem, idxglobal = BEAST.indices_splitfemglobal(gfc0 - glbf + 0.0gfcx) ## @test BEAST.Lp_integrate(gfc0 - glbf) / BEAST.Lp_integrate(glbf) ≈ 0 atol=1e-16 -@test BEAST.Lp_integrate(gfcx - glbf) / BEAST.Lp_integrate(glbf) ≈ 0.038662233651353003 +@test BEAST.Lp_integrate(gfcx - glbf) / BEAST.Lp_integrate(glbf) ≈ 0.02733832759069041 ## Global function on limited supported f(x) = x[3] diff --git a/test/test_higher_order_lagrange_functions.jl b/test/test_higher_order_lagrange_functions.jl index aec7f292..7ee6608f 100644 --- a/test/test_higher_order_lagrange_functions.jl +++ b/test/test_higher_order_lagrange_functions.jl @@ -81,8 +81,8 @@ end @test valp ≈ 1 @test crlp ≈ point(0,0,0) atol=sqrt(eps(T)) - u = T(0.2); du = eps(T) * 1000 - v = T(0.6); dv = eps(T) * 1000 + u = T(0.2); du = sqrt(eps(T)) + v = T(0.6); dv = sqrt(eps(T)) p00 = neighborhood(s, (u,v)) p10 = neighborhood(s, (u+du,v)) @@ -99,7 +99,7 @@ end for (f00, f10, f01) in zip(ϕ00, ϕ10, ϕ01) dfdu = (f10.value - f00.value)/du dfdv = (f01.value - f00.value)/dv - curl_num = (-dfdv * tu + dfdu * tv) / j + curl_num = (dfdv * tu - dfdu * tv) / j curl_ana = f00.curl @test curl_num ≈ curl_ana atol=sqrt(eps(T))*100 end diff --git a/test/test_local_assembly.jl b/test/test_local_assembly.jl index a7b5a2dd..49946454 100644 --- a/test/test_local_assembly.jl +++ b/test/test_local_assembly.jl @@ -14,7 +14,7 @@ int_nodes = submesh(!in(skeleton(boundary(m),0)), nodes) @test length(int_nodes) == 1 -L0 = BEAST.lagrangec0d1(m, int_nodes) +L0 = BEAST.lagrangec0d1(m, dirichlet=true) @test numfunctions(L0) == 1 @test length(L0.fns[1]) == 6 @test length(m) == 8 diff --git a/test/test_trace.jl b/test/test_trace.jl index f963d78c..65f69f4e 100644 --- a/test/test_trace.jl +++ b/test/test_trace.jl @@ -29,7 +29,7 @@ end @test n == 2 -## test the scalar trace of a lgrange basis +## test the scalar trace of a lagrange basis using CompScienceMeshes using BEAST using Test From 90e717e1ca3ad1dd0e11e4f241d003adf421c0d8 Mon Sep 17 00:00:00 2001 From: Cedric Muenger Date: Mon, 11 May 2026 16:32:45 +0200 Subject: [PATCH 3/6] update examples --- examples/lagc0d2.jl | 2 +- examples/mthelmholtz.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/lagc0d2.jl b/examples/lagc0d2.jl index bc39d067..ce0fc615 100644 --- a/examples/lagc0d2.jl +++ b/examples/lagc0d2.jl @@ -8,7 +8,7 @@ m = meshsphere(radius=1.0, h=0.25) nodes = skeleton(m,0) edges = skeleton(m,1) -X = BEAST.lagrangec0d2(m, nodes, edges) +X = BEAST.lagrangec0d2(m) Y = BEAST.lagrangecxd0(m) uⁱ = Helmholtz3D.planewave(wavenumber=1.0, direction=point(0,0,1)) diff --git a/examples/mthelmholtz.jl b/examples/mthelmholtz.jl index 270fbe25..491d1253 100644 --- a/examples/mthelmholtz.jl +++ b/examples/mthelmholtz.jl @@ -25,9 +25,9 @@ G12 = weld(G1,-G2) G23 = weld(G2,-G3) G31 = weld(G3,-G1) -X12 = lagrangec0d1(G12) -X23 = lagrangec0d1(G23) -X31 = lagrangec0d1(G31) +X12 = lagrangec0d1(G12; dirichlet=true) +X23 = lagrangec0d1(G23; dirichlet=true) +X31 = lagrangec0d1(G31; dirichlet=true) Y12 = duallagrangecxd0(G12) Y23 = duallagrangecxd0(G23) From 789dacd96759ae48866f32ec798cb5fd1ef5f59f Mon Sep 17 00:00:00 2001 From: Cedric Muenger Date: Tue, 12 May 2026 16:52:20 +0200 Subject: [PATCH 4/6] bugfix: interpolate type in generated lagrange polynomial expression --- src/utils/lagpolys.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/utils/lagpolys.jl b/src/utils/lagpolys.jl index e01fa77b..a6bf0406 100644 --- a/src/utils/lagpolys.jl +++ b/src/utils/lagpolys.jl @@ -35,7 +35,7 @@ end # that can be evaluated at compile time, which allows for more efficient code when the number of # nodes is known at compile time. function gen_lagpoly(nodes, i, s, i0, i1, T) - p = :(one(T)) + p = :(one($T)) for j in i0:i1 j == i && continue p = :(($p * ($s - $nodes[$j]) / ($nodes[$i] - $nodes[$j]))) @@ -44,10 +44,10 @@ function gen_lagpoly(nodes, i, s, i0, i1, T) end function gen_lagpoly_diff(nodes, i, s, i0, i1, T) - dp = :(zero(T)) + dp = :(zero($T)) for j in i0:i1 j == i && continue - p = :(one(T)) + p = :(one($T)) for k in i0:i1 (k == i || k == j) && continue p = :(($p * ($s - $nodes[$k]) / ($nodes[$i] - $nodes[$k]))) From 6d167703576c31e43163fd380357c550599d56e9 Mon Sep 17 00:00:00 2001 From: Cedric Muenger Date: Fri, 22 May 2026 13:54:45 +0200 Subject: [PATCH 5/6] remove generated shapefunctions from GWPs --- src/bases/local/gwplocal.jl | 195 ++++++++++++++++-------------------- src/utils/lagpolys.jl | 17 ++++ 2 files changed, 103 insertions(+), 109 deletions(-) diff --git a/src/bases/local/gwplocal.jl b/src/bases/local/gwplocal.jl index 296cf378..7d9f6e3f 100644 --- a/src/bases/local/gwplocal.jl +++ b/src/bases/local/gwplocal.jl @@ -22,122 +22,99 @@ function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Di ϕ(dom, u, dimtype(ϕ,dom)) end -@generated function (::GWPCurlRefSpace{T,Degree})(dom::CompScienceMeshes.ReferenceSimplex{Dim}, - bary, ::Val{NF}) where {T,Degree,Dim,NF} - - u = :(bary[1]) - v = :(bary[2]) - w = :(one(T) - $u - $v) - - #nodes for sylvester polynomials - s = :() - for i in 0:Degree+2 - s = :($s..., $i/($Degree+2)) - end - nd1 = :(SVector{2}( -$v, $u-one(T) )) - nd2 = :(SVector{2}( -$v+one(T), $u )) - nd3 = :(SVector{2}( -$v, $u )) - - nts = :() +macro gwp_shapefunction_barycentric_index(degree) + k = degree + 2 + LUT = SVector{(degree+1)*(degree+3)}(vcat( + [(0, j, k-j) for j in 1:degree+1], + [(i, 0, k-i) for i in 1:degree+1], + [(i, k-i, 0) for i in 1:degree+1], + [repeat([(i, j, k-i-j)], 2) for i in 1:degree for j in 1:degree if k-i-j ≥ 1]... + )) + esc(LUT) +end - i = 0 - for j in 1:Degree+1 - k = (Degree+2)-i-j - Rᵢ = gen_sylpoly(s, i+1, u, T) - Rⱼ = gen_sylpoly_shift(s, j+1, v, T) - Rₖ = gen_sylpoly_shift(s, k+1, w, T) - - dRᵢ = gen_sylpoly_diff(s, i+1, u, T) - dRⱼ = gen_sylpoly_shift_diff(s, j+1, v, T) - dRₖ = gen_sylpoly_shift_diff(s, k+1, w, T) - 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ₖ ) - - nts = :($nts..., (value=$Rᵢ*$Rⱼ*$Rₖ*$nd1, curl=$curl)) - end +localidx(i,::Val{0}) = @gwp_shapefunction_barycentric_index(0)[i] +localidx(i,::Val{1}) = @gwp_shapefunction_barycentric_index(1)[i] +localidx(i,::Val{2}) = @gwp_shapefunction_barycentric_index(2)[i] +localidx(i,::Val{3}) = @gwp_shapefunction_barycentric_index(3)[i] +localidx(i,::Val{4}) = @gwp_shapefunction_barycentric_index(4)[i] +localidx(i,::Val{5}) = @gwp_shapefunction_barycentric_index(5)[i] +localidx(i,::Val{6}) = @gwp_shapefunction_barycentric_index(6)[i] - for i in 1:Degree+1 - j = 0 - k = (Degree+2)-i-j - Rᵢ = gen_sylpoly_shift(s, i+1, u, T) - Rⱼ = gen_sylpoly(s, j+1, v, T) - Rₖ = gen_sylpoly_shift(s, k+1, w, T) - - dRᵢ = gen_sylpoly_shift_diff(s, i+1, u, T) - dRⱼ = gen_sylpoly_diff(s, j+1, v, T) - dRₖ = gen_sylpoly_shift_diff(s, k+1, w, T) - 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ₖ ) - - nts = :($nts..., (value=$Rᵢ*$Rⱼ*$Rₖ*$nd2, curl=$curl)) - end - for i in 1:Degree+1 - j = (Degree+2)-i - k = 0 - Rᵢ = gen_sylpoly_shift(s, i+1, u, T) - Rⱼ = gen_sylpoly_shift(s, j+1, v, T) - Rₖ = gen_sylpoly(s, k+1, w, T) - - dRᵢ = gen_sylpoly_shift_diff(s, i+1, u, T) - dRⱼ = gen_sylpoly_shift_diff(s, j+1, v, T) - dRₖ = gen_sylpoly_diff(s, k+1, w, T) - - 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ₖ ) - - nts = :($nts..., (value=$Rᵢ*$Rⱼ*$Rₖ*$nd3, curl=$curl)) - end +function shapefunction(idx::Int, bary, ::Val{Degree}) where {Degree} + + T = eltype(bary) + + u, v = bary[1], bary[2] + w = one(T)-v-u + + s = SVector{Degree+3,T}([i/(Degree+2) for i in 0:Degree+2]) - for i in 1:Degree+1 - for j in 1:Degree+1 - k = (Degree+2)-i-j - k <= 0 && continue - Rsᵢ = gen_sylpoly_shift(s, i+1, u, T) - Rsⱼ = gen_sylpoly_shift(s, j+1, v, T) - Rsₖ = gen_sylpoly_shift(s, k+1, w, T) - Rᵢ = gen_sylpoly(s, i+1, u, T) - Rⱼ = gen_sylpoly(s, j+1, v, T) - Rₖ = gen_sylpoly(s, k+1, w, T) - S1 = :( $Rᵢ*$Rsⱼ*$Rsₖ*$nd1 ) - S2 = :( $Rsᵢ*$Rⱼ*$Rsₖ*$nd2 ) - S3 = :( $Rsᵢ*$Rsⱼ*$Rₖ*$nd3 ) - - dRsᵢ = gen_sylpoly_shift_diff(s, i+1, u, T) - dRsⱼ = gen_sylpoly_shift_diff(s, j+1, v, T) - dRsₖ = gen_sylpoly_shift_diff(s, k+1, w, T) - dRᵢ = gen_sylpoly_diff(s, i+1, u, T) - dRⱼ = gen_sylpoly_diff(s, j+1, v, T) - dRₖ = gen_sylpoly_diff(s, k+1, w, T) - - 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ₖ ) - - nts = :($nts..., (value=$S2-$S3, curl=$curlS2 - $curlS3)) - - nts = :($nts..., (value=$S3-$S1, curl=$curlS3 - $curlS1)) - end end - - ex = :(SVector{NF}(())) - - for i in 1:NF - push!(ex.args[2].args, :($nts[$i])) + nd1 = SVector(-v, u-one(T)) + nd2 = SVector(-v+one(T), u) + nd3 = SVector(-v, u) + + i,j,k = localidx(idx, Val(Degree)) + + Rsᵢ = BEAST._sylpoly_shift(s, i+1, u) + Rsⱼ = BEAST._sylpoly_shift(s, j+1, v) + Rsₖ = BEAST._sylpoly_shift(s, k+1, w) + + Rᵢ = BEAST._sylpoly(s, i+1, u) + Rⱼ = BEAST._sylpoly(s, j+1, v) + Rₖ = BEAST._sylpoly(s, k+1, w) + + S1 = Rᵢ*Rsⱼ*Rsₖ*nd1 + S2 = Rsᵢ*Rⱼ*Rsₖ*nd2 + S3 = Rsᵢ*Rsⱼ*Rₖ*nd3 + + dRsᵢ = BEAST._sylpoly_shift_diff(s, i+1, u) + dRsⱼ = BEAST._sylpoly_shift_diff(s, j+1, v) + dRsₖ = BEAST._sylpoly_shift_diff(s, k+1, w) + dRᵢ = BEAST._sylpoly_diff(s, i+1, u) + dRⱼ = BEAST._sylpoly_diff(s, j+1, v) + dRₖ = BEAST._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ₖ + + if idx <= (Degree+1) + val = S1 + curl = curlS1 + elseif idx <= 2*(Degree+1) + val = S2 + curl = curlS2 + elseif idx <= 3*(Degree+1) + val = S3 + curl = curlS3 + else + if (idx-3*(Degree+1)) % 2 == 1 + val = S2 - S3 + curl = curlS2 - curlS3 + else + val = S3 - S1 + curl = curlS3 - curlS1 + end end - return ex + return (value=val, curl=curl) +end + +function (::GWPCurlRefSpace{T,Degree})(dom::CompScienceMeshes.ReferenceSimplex{Dim}, + bary, ::Val{NF}) where {T,Degree,Dim,NF} + + return SVector{NF}(shapefunction(i, bary, Val(Degree)) for i in 1:NF) end diff --git a/src/utils/lagpolys.jl b/src/utils/lagpolys.jl index a6bf0406..c3537ec6 100644 --- a/src/utils/lagpolys.jl +++ b/src/utils/lagpolys.jl @@ -31,6 +31,23 @@ function _lagpoly_diff(nodes, i, s, i0=1, i1=length(nodes)) end +function _sylpoly(nodes, i, s) + _lagpoly(nodes, i, s, 1, i) +end + +function _sylpoly_diff(nodes, i, s) + _lagpoly_diff(nodes, i, s, 1, i) +end + +function _sylpoly_shift(nodes, i, s) + _lagpoly(nodes, i, s, 2, i) +end + +function _sylpoly_shift_diff(nodes, i, s) + _lagpoly_diff(nodes, i, s, 2, i) +end + + # Versions of the above functions to be used with @generated functions. These return expressions # that can be evaluated at compile time, which allows for more efficient code when the number of # nodes is known at compile time. From 60ee1fcb3ba00fce420443bcb36252c7901344f6 Mon Sep 17 00:00:00 2001 From: Cedric Muenger Date: Fri, 22 May 2026 14:52:41 +0200 Subject: [PATCH 6/6] remove generrated shapefunction from Lagrange --- src/bases/local/gwplocal.jl | 20 ++++----- src/bases/local/laglocal.jl | 82 +++++++++++++++++-------------------- 2 files changed, 47 insertions(+), 55 deletions(-) diff --git a/src/bases/local/gwplocal.jl b/src/bases/local/gwplocal.jl index 7d9f6e3f..8044854d 100644 --- a/src/bases/local/gwplocal.jl +++ b/src/bases/local/gwplocal.jl @@ -34,16 +34,16 @@ macro gwp_shapefunction_barycentric_index(degree) esc(LUT) end -localidx(i,::Val{0}) = @gwp_shapefunction_barycentric_index(0)[i] -localidx(i,::Val{1}) = @gwp_shapefunction_barycentric_index(1)[i] -localidx(i,::Val{2}) = @gwp_shapefunction_barycentric_index(2)[i] -localidx(i,::Val{3}) = @gwp_shapefunction_barycentric_index(3)[i] -localidx(i,::Val{4}) = @gwp_shapefunction_barycentric_index(4)[i] -localidx(i,::Val{5}) = @gwp_shapefunction_barycentric_index(5)[i] -localidx(i,::Val{6}) = @gwp_shapefunction_barycentric_index(6)[i] +gwplocalidx(i,::Val{0}) = @gwp_shapefunction_barycentric_index(0)[i] +gwplocalidx(i,::Val{1}) = @gwp_shapefunction_barycentric_index(1)[i] +gwplocalidx(i,::Val{2}) = @gwp_shapefunction_barycentric_index(2)[i] +gwplocalidx(i,::Val{3}) = @gwp_shapefunction_barycentric_index(3)[i] +gwplocalidx(i,::Val{4}) = @gwp_shapefunction_barycentric_index(4)[i] +gwplocalidx(i,::Val{5}) = @gwp_shapefunction_barycentric_index(5)[i] +gwplocalidx(i,::Val{6}) = @gwp_shapefunction_barycentric_index(6)[i] -function shapefunction(idx::Int, bary, ::Val{Degree}) where {Degree} +function gwp_shapefunction(idx::Int, bary, ::Val{Degree}) where {Degree} T = eltype(bary) @@ -56,7 +56,7 @@ function shapefunction(idx::Int, bary, ::Val{Degree}) where {Degree} nd2 = SVector(-v+one(T), u) nd3 = SVector(-v, u) - i,j,k = localidx(idx, Val(Degree)) + i,j,k = gwplocalidx(idx, Val(Degree)) Rsᵢ = BEAST._sylpoly_shift(s, i+1, u) Rsⱼ = BEAST._sylpoly_shift(s, j+1, v) @@ -114,7 +114,7 @@ end function (::GWPCurlRefSpace{T,Degree})(dom::CompScienceMeshes.ReferenceSimplex{Dim}, bary, ::Val{NF}) where {T,Degree,Dim,NF} - return SVector{NF}(shapefunction(i, bary, Val(Degree)) for i in 1:NF) + return SVector{NF}(gwp_shapefunction(i, bary, Val(Degree)) for i in 1:NF) end diff --git a/src/bases/local/laglocal.jl b/src/bases/local/laglocal.jl index 5f7e5aa3..325825b3 100644 --- a/src/bases/local/laglocal.jl +++ b/src/bases/local/laglocal.jl @@ -209,63 +209,55 @@ function (f::LagrangeRefSpace{T,10,2})(mp) where T end +macro lag_shapefunction_barycentric_index(degree) + LUT = SVector{binomial(degree+2, 2)}([(i,j,k) for k in 1:degree+1 for j in 1:degree+1 for i in 1:degree+1 if i+j+k == degree+3]) + esc(LUT) +end +laglocalindex(i,::Val{0}) = @lag_shapefunction_barycentric_index(0)[i] +laglocalindex(i,::Val{1}) = @lag_shapefunction_barycentric_index(1)[i] +laglocalindex(i,::Val{2}) = @lag_shapefunction_barycentric_index(2)[i] +laglocalindex(i,::Val{3}) = @lag_shapefunction_barycentric_index(3)[i] +laglocalindex(i,::Val{4}) = @lag_shapefunction_barycentric_index(4)[i] +laglocalindex(i,::Val{5}) = @lag_shapefunction_barycentric_index(5)[i] +laglocalindex(i,::Val{6}) = @lag_shapefunction_barycentric_index(6)[i] -# Evaluate Lagrange element on a triangle with the surface curl, for arbitrary degree -@generated function (::LagrangeRefSpace{T,Degree,3,NF})(mp) where {T, Degree, NF} - - u = :(parametric(mp)[1]) - v = :(parametric(mp)[2]) - w = :(1 - $u - $v) +function lag_shapefunction(idx::Int, mp, ::Val{Degree}) where {Degree} - tu = :(tangents(mp,1)) - tv = :(tangents(mp,2)) - jd = :(jacobian(mp)) + bary = parametric(mp) + tu = tangents(mp,1) + tv = tangents(mp,2) + jd = jacobian(mp) - nodes = :() - for i in 0:Degree - nodes = :($nodes..., T($i/($Degree))) - end + T = eltype(bary) - vals = :() - diffs = :() - for k in 1:Degree+1 - for j in 1:Degree+1 - for i in 1:Degree+1 - i + j + k == Degree+3 || continue + u,v = bary[1], bary[2] + w = one(T) - u - v - pp = gen_lagpoly(nodes,i,u,1,i,T) - pq = gen_lagpoly(nodes,j,v,1,j,T) - pr = gen_lagpoly(nodes,k,w,1,k,T) - val = :($pp * $pq * $pr) + nodes = SVector{Degree+1,T}([i/(Degree) for i in 0:Degree]) - vals = :($vals..., $val) + i, j, k = laglocalindex(idx, Val(Degree)) - diffu = :(zero(T)) - diffv = :(zero(T)) + pp = _lagpoly(nodes,i,u,1,i) + pq = _lagpoly(nodes,j,v,1,j) + pr = _lagpoly(nodes,k,w,1,k) - pl = gen_lagpoly_diff(nodes,i,u,1,i,T) - diffu = :($pl*$pq*$pr) + val = pp * pq * pr - pm = gen_lagpoly_diff(nodes,j,v,1,j,T) - diffv = :($pm*$pp*$pr) - - pn = gen_lagpoly_diff(nodes,k,w,1,k,T) - diffu = :($diffu - $pn*$pp*$pq) - diffv = :($diffv - $pn*$pp*$pq) - - diffs = :($diffs..., ($diffu, $diffv)) - end - end - end + pl = _lagpoly_diff(nodes,i,u,1,i) + pm = _lagpoly_diff(nodes,j,v,1,j) + pn = _lagpoly_diff(nodes,k,w,1,k) - ex = :(SVector{NF}(())) - - for i in 1:NF - push!(ex.args[2].args, :(value=$vals[$i], curl=($diffs[$i][2]*$tu-$diffs[$i][1]*$tv)/$jd)) - end + diffu = pl*pq*pr - pn*pp*pq + diffv = pm*pp*pr - pn*pp*pq + + return (value=val, curl=(diffv*tu - diffu*tv)/jd) +end - return ex +# Evaluate Lagrange element on a triangle with the surface curl, for arbitrary degree +function (::LagrangeRefSpace{T,Degree,3,NF})(mp) where {T, Degree, NF} + + return SVector{NF}(lag_shapefunction(i, mp, Val(Degree)) for i in 1:NF) end