Skip to content

Add 1D Higher Order Lagrange Functions (Cx and C0) #175

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
krcools merged 2 commits into from
Nov 14, 2025
Merged

Add 1D Higher Order Lagrange Functions (Cx and C0) #175

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b1dc8dc
Add 1D Higher Order Lagrange Functions (Cx and C0)
sbadrian Nov 8, 2025
5ab9e12
Remove historical type instabilities
sbadrian Nov 8, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src/BEAST.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ export planewave
export RefSpace, numfunctions, coordtype, scalartype, assemblydata, geometry, refspace, valuetype
export lagrangecxd0, lagrangec0d1, duallagrangec0d1, lagrangec0d2, unitfunctioncxd0, unitfunctionc0d1
export duallagrangecxd0
export lagrangec0, lagrangecx
export lagdimension
export restrict
export raviartthomas, raowiltonglisson, positions
Expand Down
171 changes: 168 additions & 3 deletions src/bases/lagrange.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -177,7 +177,10 @@ 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 interionr 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.
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)

Expand All @@ -198,6 +201,24 @@ function lagrangec0d1_dirichlet(mesh)
lagrangec0d1(mesh, vertexlist, Val{dimension(mesh)+1})
end

function lagrangec0_dirichlet(mesh; order=1)

verts = skeleton(mesh, 0)
detached = trues(numvertices(mesh))

@krcools krcools Nov 11, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

skeleton(mesh, 0) gives the nodes (vertices that are acutally used in the cell definitions). So it is not necessarily true that numvertices(mesh) == length(skeleton(mesh,0)). The latter should be used here

@sbadrian sbadrian Nov 14, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you suggesting to replace detached = trues(numvertices(mesh)) by detached = trues(length(skeleton(mesh,0)))? I think it is correct the way it is, because the goal in the following is to filter out the detached vertices (I am following here the same recipe as at other places in lagrange.jl and I have verified that it works as desired).

for v in cells(verts)
detached[v] = false
end

bnd = boundary(mesh)
bndverts = skeleton(bnd, 0)
notonbnd = trues(numvertices(mesh))
for v in cells(bndverts)
notonbnd[v] = false
end

vertexlist = findall(notonbnd .& .!detached)

@krcools krcools Nov 11, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this can be expressed as:

setminus(skeleton(mesh,0), skeleton(boundary(mesh),0))

@sbadrian sbadrian Nov 14, 2025 •
edited
Loading

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree. In fact, the entire function could be replaced by

function lagrangec0_dirichlet(mesh; order=1)
    internalnodes = setminus(skeleton(mesh,0), skeleton(boundary(mesh),0)).cells
    internalnodes = sort!([node[1] for node in internalnodes])
    
    lagrangec0(mesh, vertexlist, Val{dimension(mesh)+1}; order=order)
end

Given that the same mechanism is used in the other lagrange functions, I could replace it everywhere in lagrange.jl.

@krcools krcools Nov 14, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's leave this all as is for now. I'll make a testitem for the case where I think the two behaviours might differ.

@krcools krcools Nov 14, 2025

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My currents view is to hide as much as possible the vertices and use the 0-skeleton instead. Sometimes meshes like gmsh add floating (=unused) vertices to the file and so the two concepts are not always interchangeable.

@sbadrian sbadrian Nov 14, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, definitely, but I think the function was written precisely in this way to filter out the unused vertices when finding the internal nodes.. I think though that the setminus approach you have proposed is indeed much more elegant and clearer. I think we just should do it consistently like this in lagrange.jl as otherwise one would wonder why there are different approaches are used.

Of course, this could be done in a separate PR as this is a general refactoring.

lagrangec0(mesh, vertexlist, Val{dimension(mesh)+1}; order=order)
end

function interior_and_junction_vertices(mesh, jct)
verts = skeleton(mesh, 0)
Expand Down Expand Up @@ -310,7 +331,8 @@ end
"""
lagrangec0d1(mesh; dirichlet=[true|false]) -> basis

Build lagrangec0d1 elements, including (dirichlet=false) or excluding (dirichlet=true) those attached to boundary vertices.
Build lagrangec0d1 elements, including (dirichlet=false)
or excluding (dirichlet=true) those attached to boundary vertices.
"""
function lagrangec0d1(mesh; dirichlet::Bool=true)
if dirichlet == false
Expand All @@ -321,6 +343,17 @@ function lagrangec0d1(mesh; dirichlet::Bool=true)
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})
Expand Down Expand Up @@ -413,8 +446,65 @@ function lagrangec0d1(mesh, vertexlist, ::Type{Val{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

function lagrangec0d1(mesh, nodes::CompScienceMeshes.AbstractMesh{U,1} where {U})
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)
Expand All @@ -438,9 +528,52 @@ function lagrangec0d1(mesh, nodes::CompScienceMeshes.AbstractMesh{U,1} where {U}
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},
Expand Down Expand Up @@ -1027,7 +1160,39 @@ 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)

Expand Down
Loading