From f27e616e6fc6e8af154ed1213a6f48e80d93e9e7 Mon Sep 17 00:00:00 2001 From: Paula Respondek Date: Mon, 28 Sep 2026 16:44:14 +0200 Subject: [PATCH 01/10] Skip parent cells with no trial cells in local refines assembly `assemble_local_refines!` maps every fine test cell to its parent cell and looks that cell up in the trial space. When the trial space is supported on only part of the parent mesh, the lookup returns 0 and `bels[0]` throws a BoundsError. `assemble_local_matched!` already guards against this. Added the same guard to `assemble_local_refines!`. Dual lagrange test functions hit this case because their support reaches into coarse cells that carry no trial functions. --- src/localop.jl | 1 + test/test_local_assembly.jl | 19 ++++++++++++++++++- 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/src/localop.jl b/src/localop.jl index ada78f1a..e2bca943 100644 --- a/src/localop.jl +++ b/src/localop.jl @@ -188,6 +188,7 @@ function assemble_local_refines!(biop::LocalOperator, tfs::Space, bfs::Space, st P = ta2g[p] Q = CompScienceMeshes.parent(tgeo, P) q = bg2a[Q] + q == 0 && continue bcell = bels[q] @assert overlap(tcell, bcell) diff --git a/test/test_local_assembly.jl b/test/test_local_assembly.jl index 49946454..866ee6e9 100644 --- a/test/test_local_assembly.jl +++ b/test/test_local_assembly.jl @@ -1,6 +1,7 @@ @info "Executing test_local_assembly.jl" using CompScienceMeshes using BEAST +using LinearAlgebra using Test @@ -46,6 +47,22 @@ BEAST.assemble_local_mixed!(Id, BC, RT, st2) Q2 = fr2() @test isapprox(Q1, Q2, atol=1e-8) +# Dual functions reach into cells of the coarse mesh that carry no trial function. +# Refines algorithm has to skip those cells. +Xd = duallagrangec0d1(m) +Yp = BEAST.subset(BEAST.lagrangecxd0(m), collect(1:2:length(m))) + +fr1, st1 = BEAST.allocatestorage(Id, Xd, Yp, Val{:bandedstorage}, BEAST.LongDelays{:ignore}) +BEAST.assemble_local_refines!(Id, Xd, Yp, st1) +Q1 = fr1() + +fr2, st2 = BEAST.allocatestorage(Id, Xd, Yp, Val{:bandedstorage}, BEAST.LongDelays{:ignore}) +BEAST.assemble_local_mixed!(Id, Xd, Yp, st2) +Q2 = fr2() + +@test isapprox(Q1, Q2, atol=1e-8) +@test norm(Q1) > 0 +@test isapprox(Matrix(assemble(Id, Xd, Yp)), Q2, atol=1e-8) @testitem "localop assembly returns sparse" begin using CompScienceMeshes @@ -71,4 +88,4 @@ Q2 = fr2() # T = typeof(Z1.lmap.A.lmap) @test Z1.lmap.A.lmap isa SparseMatrixCSC -end \ No newline at end of file +end From e193b2d7fbaea5613c2b377c0f8f5dc8a2e1fe5e Mon Sep 17 00:00:00 2001 From: Paula Respondek Date: Thu, 1 Oct 2026 16:43:38 +0200 Subject: [PATCH 02/10] Build the overlap simplex from all vertices of the child chart The rule hardcoded three vertices, so it only worked for triangles. A refined curve has two-vertex charts and raised a BoundsError. Map over the chart's vertices instead. For a triangle this produces the same three arguments as before. --- src/quadrature/rules/testrefinestrialqrule.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/quadrature/rules/testrefinestrialqrule.jl b/src/quadrature/rules/testrefinestrialqrule.jl index 0da2f922..e57f5f78 100644 --- a/src/quadrature/rules/testrefinestrialqrule.jl +++ b/src/quadrature/rules/testrefinestrialqrule.jl @@ -37,10 +37,7 @@ function integrate!(out, op, trial_charts = [chart(test_mesh, p) for p in CompScienceMeshes.children(parent_mesh, trial_cell)] trial_overlaps = map(trial_charts) do chart - simplex( - carttobary(trial_chart, chart.vertices[1]), - carttobary(trial_chart, chart.vertices[2]), - carttobary(trial_chart, chart.vertices[3])) + simplex(map(v -> carttobary(trial_chart, v), chart.vertices)) end quadstrat = qr.conforming_qstrat From bfa8805d5bf786ac46b6c6dbc2028a52da5c4b89 Mon Sep 17 00:00:00 2001 From: Paula Respondek Date: Mon, 5 Oct 2026 11:03:21 +0200 Subject: [PATCH 03/10] Remove the broken restriction of quadratic Lagrange elements The three-argument restrict for degree two called numfunctions without a domain, for which there is no method, so it threw on every call. Delete it and let the generic RefSpace fallback route through interpolate, which gives the same result as the in-place five-argument form. Add a test comparing the three routes to a restriction matrix on a genuine restriction onto a child cell. --- src/bases/local/laglocal.jl | 47 --------------------------- test/test_interpolate_and_restrict.jl | 28 ++++++++++++++++ 2 files changed, 28 insertions(+), 47 deletions(-) diff --git a/src/bases/local/laglocal.jl b/src/bases/local/laglocal.jl index 61eacffa..e9eefc60 100644 --- a/src/bases/local/laglocal.jl +++ b/src/bases/local/laglocal.jl @@ -596,53 +596,6 @@ function restrict(f::LagrangeRefSpace{T,1}, dom1, dom2) where T return Q end - -function restrict(f::LagrangeRefSpace{T,2}, dom1, dom2) where T - - D = numfunctions(f) - Q = zeros(T, D, D) - - # for each point of the new domain - for i in 1:3 - - #vertices - v = dom2.vertices[i] - - # find the barycentric coordinates in dom1 - uvn = carttobary(dom1, v) - - # evaluate the shape functions in this point - x = neighborhood(dom1, uvn) - fx = f(x) - - for j in 1:D - Q[j,i] = fx[j][1] - end - - - #edges - # find the center of edge i of dom2 - a = dom2.vertices[mod1(i+1,3)] - b = dom2.vertices[mod1(i+2,3)] - v = (a + b) / 2 - - # find the barycentric coordinates in dom1 - uvn = carttobary(dom1, v) - - # evaluate the shape functions in this point - x = neighborhood(dom1, uvn) - fx = f(x) - - for j in 4:D - Q[j,i+3] = fx[j][1] - end - end - - return Q -end - - - # Quadratic Lagrange element on a triangle # function (f::LagrangeRefSpace{T,2,3})(t) where T # u,v,w, = barycentric(t) diff --git a/test/test_interpolate_and_restrict.jl b/test/test_interpolate_and_restrict.jl index 9671c0e6..255fbafc 100644 --- a/test/test_interpolate_and_restrict.jl +++ b/test/test_interpolate_and_restrict.jl @@ -58,4 +58,32 @@ itpol = sum(w*val for (w,val) in zip(Q3,vals)) end end +@testitem "restrict Lagrange on a refined cell" begin + using CompScienceMeshes + using LinearAlgebra + + coarse = simplex(point(0,0,0), point(1,0,0), point(0,1,0)) + # one child of the barycentric refinement of coarse: + child = simplex(point(1,0,0),point(0.5,0.5,0), point(1/3,1/3,0)) + + tobary(chart, ch) = simplex(map(v -> carttobary(chart,v), ch.vertices)) + # the three routes to a restriction matrix must agree on a genuine + # restriciton; self_restriction alone passes under any node ordering + for (degree, n) in ((1,3),(2,6)) + rs = BEAST.LagrangeRefSpace{Float64, degree, 3, n}() + Q1 = BEAST.restrict(rs, coarse, child) + Q2 = BEAST.restrict(rs, coarse, child, tobary(coarse, child)) + Q3 = zeros(n, n) + BEAST.restrict!(Q3, rs, coarse, child, tobary(coarse, child)) + @test Q1 ≈ Q2 + @test Q1 ≈ Q3 + end + + for (degree, n) in ((0,1),(1,3),(2,6),(3,10)) + rs = BEAST.LagrangeRefSpace{Float64, degree, 3, n}() + Q = zeros(n, n) + BEAST.restrict!(Q, rs, coarse, coarse, tobary(coarse, coarse)) + @test Q ≈ I + end +end From cf256c32bd6658562e1747fd303f2f00dde7a0f8 Mon Sep 17 00:00:00 2001 From: Paula Respondek Date: Wed, 30 Sep 2026 14:13:14 +0200 Subject: [PATCH 04/10] Fix duallagrangecxd0 on curves For a curve the default junction is zero-dimensional, so the call dispatched to the method that reads its second argument as a vertex list given as a mesh. The empty default junction became an empty vertex list and the constructor silently returned a space with no functions. A junction passed explicitly hit the same method, yielding dual functions on the two end vertices only instead of on every vertex. Give the one-argument form its own method and decide between junction and vertex list by comparing dimensions rather than by dispatch, which cannot separate the two readings on a curve. Surfaces are unaffected. --- src/bases/lagrange.jl | 46 ++++++++++++++++++++++++++++++++----------- 1 file changed, 34 insertions(+), 12 deletions(-) diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl index d90abca9..58fe0a4d 100644 --- a/src/bases/lagrange.jl +++ b/src/bases/lagrange.jl @@ -814,11 +814,20 @@ the area so that overall the integral over the dual function is one. When `interpolatory=true` is used, the function value is one on the support, and thus, it gives rise to a partition of unity. """ -function duallagrangecxd0(mesh, jct=CompScienceMeshes.mesh(coordtype(mesh), dimension(mesh)-1); interpolatory=false) - vertexlist = interior_and_junction_vertices(mesh, jct) - duallagrangecxd0(mesh, vertexlist; interpolatory=interpolatory) +function duallagrangecxd0(mesh; interpolatory=false) + jct = CompScienceMeshes.mesh(coordtype(mesh), dimension(mesh)-1) + duallagrangecxd0(mesh, jct; interpolatory=interpolatory) end +function duallagrangecxd0(mesh, jct::CompScienceMeshes.AbstractMesh; interpolatory=false) + if dimension(jct) == dimension(mesh) - 1 + vertexlist = interior_and_junction_vertices(mesh, jct) + else + @assert dimension(jct) == 0 + vertexlist = Int[CompScienceMeshes.indices(jct, v)[1] for v in jct] + end + duallagrangecxd0(mesh, vertexlist; interpolatory=interpolatory) +end function duallagrangecxd0(mesh, vertexlist::Vector{Int}; interpolatory=false) @@ -841,14 +850,6 @@ function duallagrangecxd0(mesh, vertexlist::Vector{Int}; interpolatory=false) LagrangeBasis{0,-1,NF}(fine, fns, pos) end - -function duallagrangecxd0(mesh, vertices::CompScienceMeshes.AbstractMesh{U,1}; interpolatory=false) where {U} - # vertexlist = Int[v[1] for v in vertices] - vertexlist =Int[CompScienceMeshes.indices(vertices, v)[1] for v in vertices] - return duallagrangecxd0(mesh, vertexlist; interpolatory=interpolatory) -end - - """ singleduallagd0(fine, F, v; interpolatory=false) @@ -1962,4 +1963,25 @@ end using CompScienceMeshes Γ = meshcuboid(1.0,1.0,1.0,0.5) @test BEAST._surface(duallagrangecxd0(Γ)) ≈ 6.0 -end \ No newline at end of file +end + +@testitem "duallagrangecxd0 on curves" begin + using CompScienceMeshes + + # closed curve, vertex numbering inherited from the parent surface + Γ = boundary(meshrectangle(1.0,1.0,0.5,3)) + X = duallagrangecxd0(Γ) + @test numfunctions(X) == numcells(Γ) + @test BEAST._surface(X) ≈ 4.0 + + # open curve: the two end vertices carry no dual function + Λ = meshsegment(1.0, 1/4, 3) + Y = duallagrangecxd0(Λ) + @test numfunctions(Y) == numcells(Λ) - 1 + @test BEAST._surface(Y) ≈ 0.75 + + # junction supplied, half-pyramids on end vertices are kept + Z = duallagrangecxd0(Λ, boundary(Λ)) + @test numfunctions(Z) == numcells(Λ) + 1 + @test BEAST._surface(Z) ≈ 1.0 +end From 6ee614313e6625fcc6f5bcf3c326ae652a26cb1a Mon Sep 17 00:00:00 2001 From: Paula Respondek Date: Thu, 1 Oct 2026 14:25:09 +0200 Subject: [PATCH 05/10] Generalize the dual linear Lagrange coefficients on surfaces The three coefficients were hardcoded for a closed surface: one at the barycentre, one half at a face centre, and one over half the number of fine cells at a coarse vertex. The face-centre value assumes two coarse cells share every edge, which fails on an open surface, where a boundary edge belongs to one. The dual functions there summed to one half rather than one, so the space could not represent a constant. Replace all three with the fraction of the fine cells meeting at a vertex that belong to the coarse cell, and fold the three loops into one. On a closed surface this reproduces the previous shapes exactly. --- src/bases/lagrange.jl | 51 +++++++++++++++++-------------------------- 1 file changed, 20 insertions(+), 31 deletions(-) diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl index 58fe0a4d..61ce0622 100644 --- a/src/bases/lagrange.jl +++ b/src/bases/lagrange.jl @@ -1193,6 +1193,7 @@ function duallagrangec0d1(mesh, refined, jct_pred, ::Type{Val{3}}) uv_ctr = ones(dimension(mesh))/(dimension(mesh)+1) vtoc, vton = vertextocellmap(refined) + fine_cells = collect(refined) for (i,p) in enumerate(mesh) coarse_idcs = CompScienceMeshes.indices(mesh, p) coarse_chart = chart(mesh,p) @@ -1209,6 +1210,7 @@ function duallagrangec0d1(mesh, refined, jct_pred, ::Type{Val{3}}) I = CollisionDetection.find(fine_vertices, centroid) @assert length(I) == 1 centroid_id = I[1] + children = vtoc[centroid_id, 1:vton[centroid_id]] # get the indx in fine.vertices of the centroid of the faces of coarse_cell face_center_ids = Vector{Int}(undef,num_faces) @@ -1225,39 +1227,16 @@ function duallagrangec0d1(mesh, refined, jct_pred, ::Type{Val{3}}) face_center_ids[f] = I[1] end - n = vton[centroid_id] - for c in vtoc[centroid_id,1:n] - fine_idcs = cells(refined)[c].indices - local_id = something(findfirst(isequal(centroid_id), fine_idcs), 0) - @assert local_id != 0 - shape = Shape(c, local_id, 1.0) - push!(fns[i], shape) - end - - for f in 1:num_faces - v = face_center_ids[f] - jct_pred(vertices(refined)[v]) && continue - n = vton[v] - for c in vtoc[v,1:n] - fine_idcs = cells(refined)[c].indices - local_id = something(findfirst(isequal(v), fine_idcs),0) - @assert local_id != 0 - #shape = Shape(c, local_id, 1/n/2) - shape = Shape(c, local_id, 1/2) - - push!(fns[i], shape) - end - end - - for f in 1:length(coarse_idcs) - v = coarse_idcs[f] - jct_pred(vertices(refined)[v]) && continue + for v in (centroid_id, face_center_ids..., coarse_idcs...) + v != centroid_id && jct_pred(vertices(refined)[v]) && continue n = vton[v] - for c in vtoc[v,1:n] - fine_idcs = cells(refined)[c].indices - local_id = something(findfirst(isequal(v), fine_idcs),0) + cells_at_v = vtoc[v, 1:n] + coeff = T(count(in(children), cells_at_v)) / n + for c in cells_at_v + fine_idcs = CompScienceMeshes.indices(refined, fine_cells[c]) + local_id = something(findfirst(isequal(v), fine_idcs), 0) @assert local_id != 0 - shape = Shape(c, local_id, 1/(n/2)) + shape = Shape(c, local_id, coeff) push!(fns[i], shape) end end @@ -1959,6 +1938,16 @@ end @test BEAST._surface(duallagrangec0d1(Γ)) ≈ 6.0 end +@testitem "duallagrangec0d1 on open surfaces" begin + using CompScienceMeshes + # dual P1 must reproduce constants on an open surface too: a boundary edge + # belongs to one coarse cell, so its coefficient is one, not one half + Γ = meshrectangle(1.0, 1.0, 0.25, 3) + @test BEAST._surface(duallagrangec0d1(Γ)) ≈ 1.0 + # closed surfaces are unchanged + @test BEAST._surface(duallagrangec0d1(meshcuboid(1.0,1.0,1.0,0.5))) ≈ 6.0 +end + @testitem "duallagrangecxd0" begin using CompScienceMeshes Γ = meshcuboid(1.0,1.0,1.0,0.5) From 3106d7c273211860ca4676528e2564cdd28274bd Mon Sep 17 00:00:00 2001 From: Paula Respondek Date: Mon, 5 Oct 2026 12:39:30 +0200 Subject: [PATCH 06/10] Rewrite dual linear Lagrange elements on curves Upstream fixed the array size, but the shape assignment itself was wrong: The construction assigned the shape at a coarse vertex only on the fine cell inside the coarse segment, and put a shape at the neighbouring segment's midpoint on the adjacent cell. The resulting functions jumped by one across every coarse vertex, so the space was not continuous. It also indexed of a surface, indexed a second adjacent cell that does not exist at the end of an open curve, reached into mesh.faces directly, and ignored the junction predicate it was passed. Assign to every fine vertex the fraction of the fine cells meeting there that belong to the coarse segment, on each incident cell, as the surface construction above now does. --- src/bases/lagrange.jl | 113 ++++++++++++++++-------------------------- 1 file changed, 42 insertions(+), 71 deletions(-) diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl index 61ce0622..e7d2b371 100644 --- a/src/bases/lagrange.jl +++ b/src/bases/lagrange.jl @@ -1259,79 +1259,50 @@ This basis function creats the dual Lagrange basis function and return an object It also return a gemoetry containing the refined mesh. """ function duallagrangec0d1(mesh, mesh2, pred, ::Type{Val{2}}) - T = coordtype(mesh) - U = universedimension(mesh) - # get the information about number of vertices, number of faces , and the maping between vertices and faces for the original mesh - numverts1 = numvertices(mesh) - num_cells1 = numcells(mesh) - cellids1, ncells1=vertextocellmap(mesh) - # get the information about number of vertices, number of faces , and the maping between vertices and faces for the refined mesh - num_cells2 = numcells(mesh2) - numverts2 = numvertices(mesh2) - geometry = mesh2 - cellids2, ncells2 = vertextocellmap(mesh2) - - fns = Vector{Vector{Shape{T}}}() - pos = Vector{vertextype(mesh)}() - # We will iterate over the coarse mesh segments to assign all the functions to it. - for segment_coarse in 1 : num_cells1 - # For the dual Lagrange there is a 6 shapes per segment - numshapes = (ncells1[mesh.faces[segment_coarse].indices[1]]*4) -2 - shapes = Vector{Shape{T}}(undef,numshapes) - # Now we will get all the smaller faces within the coarse segment - #i.e The coose segment will have two points, and these tow points are connected to two segments in the finer mesh - # This will give us a 4 smaller faces per Dual lagrange basis, we store them first in all_faces - # all_faces= Array{SVector{2,Int}}(undef,4) # faces in the original segment (4) - CT = CompScienceMeshes.SimplexGraph{2} - all_faces = Vector{CT}(undef, 4) - # the follwoing code get the verteciec for each coarse segment - # then it looks for the two faces connected to each point in the finer mesh - # if for example the segment would connect to more than two faces we will have - # mesh.faces[segment_coarse][1],n] and iterate over how many segments [n] are connected - all_faces[1]=mesh2.faces[cellids2[mesh.faces[segment_coarse][1],1]] - all_faces[2]=mesh2.faces[cellids2[mesh.faces[segment_coarse][1],2]] - all_faces[3]=mesh2.faces[cellids2[mesh.faces[segment_coarse][2],1]] - all_faces[4]=mesh2.faces[cellids2[mesh.faces[segment_coarse][2],2]] - - # now we now the first point of the corse segment will have the left hand side basis - # and we know that that point is connected to faces 1,2 in the array all_faces - # but now we want to which one of them is inner segment and which one is at the edge - # as both of these two faces are in left hand side - # The inner face will have the left hand side corase point in its first place for example - # we have original point 2,3 , and the coarse segment (2,3) - # after refinment we will have for example points 20,2,21,3,22 with four faces (20,2),(2,21),(21,3),(3,22) - # so the faces to the left (20,2),(2,21) we can decide the inner one if the original point(2) is at first index - # same for right hand side faces (21,3),(3,22) if the original point(3) is in the second index - # then the inner faces are (2,21) and (21,3) - # the following code does the same checking process and assign the shapes for the dual lagragne right away - - # For the Left hand side faces - if(all_faces[1][1] == mesh.faces[segment_coarse][1]) - shapes[1]= Shape(cellids2[mesh.faces[segment_coarse][1],2],1,0.5) - shapes[2]= Shape(cellids2[mesh.faces[segment_coarse][1],1],2,0.5) - shapes[3]= Shape(cellids2[mesh.faces[segment_coarse][1],1],1,1.0) - elseif(all_faces[2][1] == mesh.faces[segment_coarse][1]) - shapes[1]= Shape(cellids2[mesh.faces[segment_coarse][1],2],1,1.0) - shapes[2]= Shape(cellids2[mesh.faces[segment_coarse][1],2],2,0.5) - shapes[3]= Shape(cellids2[mesh.faces[segment_coarse][1],1],1,0.5) - end - # For the Right hand side faces - if(all_faces[3][2] == mesh.faces[segment_coarse][2]) - shapes[4]= Shape(cellids2[mesh.faces[segment_coarse][2],1],2,1.0) - shapes[5]= Shape(cellids2[mesh.faces[segment_coarse][2],1],1,0.5) - shapes[6]= Shape(cellids2[mesh.faces[segment_coarse][2],2],2,0.5) - elseif(all_faces[4][2] == mesh.faces[segment_coarse][2]) - shapes[4]= Shape(cellids2[mesh.faces[segment_coarse][2],2],2,1.0) - shapes[5]= Shape(cellids2[mesh.faces[segment_coarse][2],2],1,0.5) - shapes[6]= Shape(cellids2[mesh.faces[segment_coarse][2],1],2,0.5) + T = coordtype(mesh) + + fns = Vector{Vector{Shape{T}}}(undef, numcells(mesh)) + pos = Vector{vertextype(mesh)}() + + # store the fine mesh's vertices in an octree for fast retrieval + fine_vertices = Octree(vertices(mesh2)) + vtoc, vton = vertextocellmap(mesh2) + fine_cells = collect(mesh2) + + for (i, p) in enumerate(mesh) + coarse_idcs = CompScienceMeshes.indices(mesh, p) + centroid = cartesian(center(chart(mesh, p))) + + fns[i] = Vector{Shape{T}}() + push!(pos, centroid) + + # midpoint of a coarse segment is a vertex of the refinement, and the + # fine cells meeting there are the two children of that segment + I = CollisionDetection.find(fine_vertices, centroid) + @assert length(I) == 1 + centroid_id = I[1] + children = vtoc[centroid_id, 1:vton[centroid_id]] + + # The coefficient at a fine vertex is the fraction of the fine cells meeting + # there that belong to this coarse segment: one at the midpoint, one half at + # an interior coarse vertex, one at the end of an open curve. + for v in (centroid_id, coarse_idcs...) + v != centroid_id && pred(vertices(mesh2)[v]) && continue + n = vton[v] + cells_at_v = vtoc[v, 1:n] + coeff = T(count(in(children), cells_at_v)) / n + for c in cells_at_v + fine_idcs = CompScienceMeshes.indices(mesh2, fine_cells[c]) + local_id = something(findfirst(isequal(v), fine_idcs), 0) + @assert local_id != 0 + push!(fns[i], Shape(c, local_id, coeff)) + end + end end - # Now assign all of these shapes to the relevent segment in the coarse mesh - push!(fns, shapes) - push!(pos, cartesian(center(chart(mesh, segment_coarse)))) - end + @assert length(fns) == length(pos) - NF = 2 - LagrangeBasis{1,0,NF}(geometry, fns, pos) + NF = 2 + LagrangeBasis{1,0,NF}(mesh2, fns, pos) end # """ From ac51bfed1286970c23cc755f70d0c5b7da788ea7 Mon Sep 17 00:00:00 2001 From: Paula Respondek Date: Mon, 5 Oct 2026 13:21:20 +0200 Subject: [PATCH 07/10] Add interpolation for segments and tetrahedra Lagrange reference spaces on segments and tetrahedra had no interpolat!, so restrict! thre a MethodError for every one-dimensional cell and tetrahedral degreee one. Degree zero on a tetrahedron had no evaluator at all. Add the missing methods in the same inline style as the triangle versions. The segment nodes are the two end points followed by the interior nodes in descending order, which is the order the hand-written evaluators one to ten use. On a tetrahedron degree one interpolates at the four vertices; higher degrees have no evaluator, so nothing further is needed at the moment. --- src/bases/local/laglocal.jl | 59 ++++++++++++++++++++++++++- test/test_interpolate_and_restrict.jl | 27 ++++++++++++ 2 files changed, 84 insertions(+), 2 deletions(-) diff --git a/src/bases/local/laglocal.jl b/src/bases/local/laglocal.jl index e9eefc60..eca8e224 100644 --- a/src/bases/local/laglocal.jl +++ b/src/bases/local/laglocal.jl @@ -15,7 +15,7 @@ valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} = T # Evaluate constant lagrange elements on anything (ϕ::LagrangeRefSpace{T,0,2})(tp) where {T} = SVector(((value=one(T), derivative=zero(T)),)) - +(ϕ:: LagrangeRefSpace{T,0,4})(tp) where {T} = SVector(((value=one(T), derivative=zero(T)),)) # Evaluate linear Lagrange elements on a segment # The derivative denotes the tangential derivative function (f::LagrangeRefSpace{T,1,2})(mp) where T @@ -858,4 +858,59 @@ function interpolate!(out, fields, interpolant::LagrangeRefSpace{T,0,3}, chart) vals = fields(p) for (g, val) in zip(axes(out, 1), vals) out[g,1] = val -end end \ No newline at end of file +end end + +# segment methods +function interpolate(fields, interpolant::LagrangeRefSpace{T, Degree, 2},chart) where {T, Degree} + vals = Vector{Vector{T}}() + if Degree > 0 + push!(vals, fields(neighborhood(chart, (T(1),)))) + push!(vals, fields(neighborhood(chart, (T(0),)))) + for k in 1:Degree-1 + push!(vals, fields(neighborhood(chart, (T(Degree-k)/Degree,)))) + end + else + push!(vals, fields(center(chart))) + end + 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, 2}, chart) where {T, Degree} + us = Degree > 0 ? (T(1), T(0), (T(Degree-k)/Degree for k in 1:Degree-1)...) : (T(1)/2,) + for (idx, u) in enumerate(us) + vals = fields(neighborhood(chart, (u,))) + for (g, val) in zip(axes(out, 1), vals) + out[g, idx] = val + end + end +end + +# tetrahedron methods +function interpolate(fields, interpolant::LagrangeRefSpace{T, 1, 4}) where {T} + vals = [fields(neighborhood(chart, uvw)) for uvw in ((T(1),T(0),T(0)),(T(0),T(1),T(0)),(T(0),T(0),T(1)),(T(0),T(0),T(0)))] + 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, 1, 4}, chart) where {T} + for (idx, uvw) in enumerate(((T(1),T(0),T(0)),(T(0),T(1),T(0)),(T(0),T(0),T(1)),(T(0),T(0),T(0)))) + vals = fields(neighborhood(chart, uvw)) + for (g, val) in zip(axes(out, 1), vals) + out[g, idx] = val + end + end +end + +function interpolate!(out, fields, interpolant::LagrangeRefSpace{T,0,4}, chart) where {T} + vals = fields(center(chart)) + for (g, val) in zip(axes(out, 1), vals) + out[g, 1] = val + end +end diff --git a/test/test_interpolate_and_restrict.jl b/test/test_interpolate_and_restrict.jl index 255fbafc..bfa53bde 100644 --- a/test/test_interpolate_and_restrict.jl +++ b/test/test_interpolate_and_restrict.jl @@ -86,4 +86,31 @@ end BEAST.restrict!(Q, rs, coarse, coarse, tobary(coarse, coarse)) @test Q ≈ I end + + # segments: nodes are the two ends followed by interior points descending + seg = simplex(point(1,0,0), point(0,0,0)) + for (degree, n) in ((0,1),(1,2),(2,3),(3,4),(4,5)) + rs = BEAST.LagrangeRefSpace{Float64, degree, 2, n}() + Q = zeros(n, n) + BEAST.restrict!(Q, rs, seg, seg, tobary(seg, seg)) + @test Q ≈ I + end + + # genuine restriction on a segment: child is the left half + half = simplex(point(1,0,0), point(0.5,0,0)) + for ( degree, n) in ((1,2),(2,3),(3,4)) + rs = BEAST.LagrangeRefSpace{Float64, degree, 2, n}() + Q = zeros(n, n) + BEAST.restrict!(Q, rs, seg, half, tobary(seg, half)) + @test Q ≈ BEAST.restrict(rs, seg, half) + end + + # tetrahedra: only degrees zero and one have an evaluator + tet = simplex(point(0,0,0), point(1,0,0),point(0,1,0), point(0,0,1)) + for (degree, n) in ((0,1),(1,4)) + rs = BEAST.LagrangeRefSpace{Float64, degree, 4, n}() + Q = zeros(n, n) + BEAST.restrict!(Q, rs, tet, tet, tobary(tet, tet)) + @test Q ≈ I + end end From 48f5064cec1d59301c635f3d153f66b4bc4be837 Mon Sep 17 00:00:00 2001 From: Paula Respondek Date: Fri, 2 Oct 2026 11:22:28 +0200 Subject: [PATCH 08/10] Fix unitfunctioncxd0 in a two dimensional universe The average position was built from three literal coordinates, so every mesh whose vertices have two components raised a DimensionMismatch. Use zero(vertextype(mesh)). Makes _surface usable for curves and surfaces in the plane. Also add the testitem for dual linear Lagrange on curves. It covers a closed polygon, an open segment, a curve whose vertex numbering is inherited from a parent surface, and the junction variant, asserting continuity and the shape counts at open ends. --- src/bases/lagrange.jl | 66 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 65 insertions(+), 1 deletion(-) diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl index e7d2b371..81e9e4ec 100644 --- a/src/bases/lagrange.jl +++ b/src/bases/lagrange.jl @@ -601,7 +601,7 @@ function unitfunctioncxd0(mesh) # in case of a global function. Might be replaced by something # more useful. # For now, we fill it with the average position of the shape functions - p = vertextype(mesh)(0.0, 0.0, 0.0) + p = zero(vertextype(mesh)) for cell in mesh p += cartesian(center(chart(mesh, cell))) end @@ -1945,3 +1945,67 @@ end @test numfunctions(Z) == numcells(Λ) + 1 @test BEAST._surface(Z) ≈ 1.0 end + +@testitem "duallagrangec0d1 on curves" begin + using CompScienceMeshes + using StaticArrays + + function maxjump(X) #function to find maximum jump of basis X + fine = BEAST.geometry(X) + fcells = collect(fine) + incident = Dict{Int, Vector{Int}}() + for c in 1:numcells(fine), v in CompScienceMeshes.indices(fine, fcells[c]) + push!(get!(incident, v, Int[]), c) + end + jump = 0.0 + for fn in X.fns + val = Dict{Tuple{Int,Int},Float64}() + for s in fn + v = CompScienceMeshes.indices(fine, fcells[s.cellid])[s.refid] + val[(s.cellid, v)] = get(val, (s.cellid, v), 0.0) + s.coeff + end + for (c, v) in keys(val) + xs = [get(val, (c2, v), 0.0) for c2 in incident[v]] + jump = max(jump, maximum(xs) - minimum(xs)) + end + end + jump + end + + len(m) = sum(volume(chart(m, p)) for p in m) + + #close curve + n = 8 + verts = [SVector(cos(2π*k/n), sin(2π*k/n), 0.0) for k in 0:n-1] + faces = [CompScienceMeshes.SimplexGraph{2}(SVector(k, mod1(k+1, n))) for k in 1:n] + Γ = Mesh(verts, faces) + X = duallagrangec0d1(Γ) + @test numfunctions(X) == numcells(Γ) + @test all(==(6), length.(X.fns)) + @test maxjump(X) < 1e-12 # check if space is C0 + @test BEAST._surface(X) ≈ len(Γ) + + # open curve: end segments carry one shape less + Λ = meshsegment(1.0,1/4, 3) + Y = duallagrangec0d1(Λ) + @test numfunctions(Y) == numcells(Λ) + @test sort(unique(length.(Y.fns))) == [5, 6] + @test maxjump(Y) < 1e-12 + @test BEAST._surface(Y) ≈ 1.0 + + # vertex numbering inherited from a parent surface (github issue #195) + Γb = boundary(meshrectangle(1.0,1.0,1/4, 2)) + Z = duallagrangec0d1(Γb) + @test numfunctions(Z) == numcells(Γb) + @test maxjump(Z) < 1e-12 +end + +@testitem "unitfunctioncxd0 in a two dimensional universe" begin + using CompScienceMeshes + + Γ = meshcircle(1.0, 2π/20) # a curve in the plane + len = sum(volume(chart(Γ, p)) for p in Γ) + @test numfunctions(BEAST.unitfunctioncxd0(Γ)) == 1 + @test BEAST._surface(duallagrangec0d1(Γ)) ≈ len + @test BEAST._surface(duallagrangecxd0(Γ)) ≈ len +end From d2a642ebef2316ed9dcd2f194c5f96319f07ed0a Mon Sep 17 00:00:00 2001 From: Paula Respondek Date: Fri, 2 Oct 2026 10:29:27 +0200 Subject: [PATCH 09/10] Test refinement aware quadrature on a curve Pair a dual linear Lagrange space on the barycentric refinement of a circle with piecewise constants on the coarse curve, and compare the assembled single layer against the same operator with both spaces on the fine mesh. Before the refinement of a curve was tagged, this pairing fell back to conforming quadrature and produced 63 infinite entries out of 441. The test needs the CompScienceMeshes change that returns a BarycentricRefinement for curves, so it fails against released version (currently 0.12.1) until that is merged. --- src/quadrature/rules/testrefinestrialqrule.jl | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/quadrature/rules/testrefinestrialqrule.jl b/src/quadrature/rules/testrefinestrialqrule.jl index e57f5f78..7735c02e 100644 --- a/src/quadrature/rules/testrefinestrialqrule.jl +++ b/src/quadrature/rules/testrefinestrialqrule.jl @@ -61,3 +61,29 @@ function integrate!(out, op, for k in 1:size(Q, 2) out[i,j] += zlocal[i,k] * Q[j,k] end end end end end + +@testitem "TestRefinesTrialQRule on a curve" begin + using BEAST, Test + using CompScienceMeshes + using LinearAlgebra + + Γ = meshcircle(1.0, 2π/20) + X = duallagrangec0d1(Γ) + Y = lagrangecxd0(Γ) + @test CompScienceMeshes.refines(BEAST.geometry(X), Γ) + + 𝒱 = Helmholtz2D.singlelayer(wavenumber=1.0) + A = Matrix(assemble(𝒱, X, Y)) + @test all(isfinite, A) + + # reference: lift the coarse constants onto the fine mesh, so both spaces live on the + # same geometry and the conforming rules apply + fine = BEAST.geometry(X) + fns = [[BEAST.Shape(2*(i-1)+1, 1, 1.0), BEAST.Shape(2*(i-1)+2, 1, 1.0)] for i in 1:numcells(Γ)] + pos = [cartesian(CompScienceMeshes.center(chart(Γ,p))) for p in Γ] + Yf = BEAST.LagrangeBasis{0,-1,1}(fine, fns, pos) + @test Matrix(assemble(BEAST.Identity(), X, Yf)) ≈ Matrix(assemble(BEAST.Identity(), X, Y)) + + B = Matrix(assemble(𝒱, X, Yf)) + @test norm(A-B)/norm(B) < 1e-12 +end From 7c67043867d6df740ef1195886ef1a8600525ecc Mon Sep 17 00:00:00 2001 From: Paula Respondek Date: Fri, 2 Oct 2026 11:25:48 +0200 Subject: [PATCH 10/10] Add test for TrialRefinesTestQRule on curve --- src/quadrature/rules/testrefinestrialqrule.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/quadrature/rules/testrefinestrialqrule.jl b/src/quadrature/rules/testrefinestrialqrule.jl index 7735c02e..35d36856 100644 --- a/src/quadrature/rules/testrefinestrialqrule.jl +++ b/src/quadrature/rules/testrefinestrialqrule.jl @@ -86,4 +86,10 @@ end end end end end B = Matrix(assemble(𝒱, X, Yf)) @test norm(A-B)/norm(B) < 1e-12 + + # opposite orientation goes through TrialRefinesTestQRule + A2 = Matrix(assemble(𝒱, Y, X)) + B2 = Matrix(assemble(𝒱, Yf, X)) + @test all(isfinite,A2) + @test norm(A2 - B2)/norm(B2) < 1e-12 end