From 93f42d99bf284dcf53fbb4f618aa5186f0a261c4 Mon Sep 17 00:00:00 2001 From: Paula Respondek Date: Thu, 1 Oct 2026 15:19:43 +0200 Subject: [PATCH 1/2] Fix isinclosure for zero-dimensional cells The type for the tolerance came from the barycentric coordinates of the point, which are empty for a cell of dimension zero, so eltype returned Union{} and eps threw a MethodError. Now take the type from the cell instead. The rest of the function already treats a point cell correctly. This makes inclosure_gpredicate usable on the boundary of a curve, which is a mesh of nodes. --- src/isinside.jl | 2 +- test/test_isinclosure.jl | 9 +++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src/isinside.jl b/src/isinside.jl index eb4bca1..9203676 100644 --- a/src/isinside.jl +++ b/src/isinside.jl @@ -29,7 +29,7 @@ Determine whether point is in the closure of simplex. False positives are possib function isinclosure(cell, point) u = carttobary(cell, point) - T = eltype(u) + T = coordtype(cell) dim = dimension(cell) udim = length(point) diff --git a/test/test_isinclosure.jl b/test/test_isinclosure.jl index 3fe0ab1..b340d87 100644 --- a/test/test_isinclosure.jl +++ b/test/test_isinclosure.jl @@ -15,3 +15,12 @@ f = inclosure_gpredicate(m) @test f(p2) == false @test f(p3) == true @test f(p4) == false + +## zero-dimensional cells: the closure of a point is the point itself +b = boundary(meshsegment(1.0,1/4,3)) +g = inclosure_gpredicate(b) + +@test g(point(0.0,0,0)) == true #start point +@test g(point(1.0,0,0)) == true # end point +@test g(point(0.5,0,0)) == false # interior point +@test g(point(0.0,1.0,0)) == false # point off the segment From 6639785ccde7d544ae0139c8fe2f981d3a2efc40 Mon Sep 17 00:00:00 2001 From: Paula Respondek Date: Fri, 2 Oct 2026 09:09:37 +0200 Subject: [PATCH 2/2] Return a BarycentricRefinement when refining a curve Refining a curve returned a plain Mesh, so refines reported false and callers pairing a dual space on the refinement with a primal space on the coarse curve silently fell back to conforming quadrature, which applies singular rules to overlapping segments. Build the refinement through the generic mesh interface so it also accepts a submesh, return a BarycentricRefinement as the surface and volume versions do, and add the parent override, since the generic one returns the child index rather than the parent cell. --- src/baryref.jl | 47 +++++++++++++++++++------------------------- test/test_baryref.jl | 29 ++++++++++++++++++++++++++- 2 files changed, 48 insertions(+), 28 deletions(-) diff --git a/src/baryref.jl b/src/baryref.jl index b218ebd..a1021b1 100644 --- a/src/baryref.jl +++ b/src/baryref.jl @@ -71,52 +71,43 @@ Create the mesh obtained by inserting an extra vertex in the barycenters of all recusively creating fine cells by connecting the barycenter of a k-cell to the already constructed refined (k-1)-cells on its boundary. """ -function barycentric_refinement(mesh::Mesh{U,2}; sort=:spacefillingcurve) where U - +function barycentric_refinement(mesh::AbstractMesh{U,2}; sort=:spacefillingcurve) where U + T = coordtype(mesh) # Get the points and the faces (segments) we will use in the refinement process Edges = skeleton(mesh, 1; sort) - # Note thier number and use it to find the new number of points and segments(faces) - # after refienmnts - NV, NE = numvertices(Edges), numcells(Edges) - - nv = NV + NE - ne = 2NE + # Note their number and use it to size the new vertex and segment arrays + NV, NE = numvertices(mesh), numcells(Edges) # define the new vertices array that will hold the old points(coarse) and the # new points(finer mesh) of vertices - verts = similar(mesh.vertices, nv) + verts = Array{vertextype(mesh)}(undef, NV+NE) - # Now we assign the coarse mesh vetrices (old points) in the upper half excatly - # as the original coarse mesh + # Now we assign the coarse mesh vetrices (old points) to the first half of the + # array excatly as the original coarse mesh for i in 1:NV - verts[i] = mesh.vertices[i] - end - - # for the second half we store the refienment points and we calculate them using - # the old points. We insert a point between every pair of points in the original - # coarse mesh by (point1+ point2)\2 - # for (i, Face) in enumerate(cells(mesh)) - for (i,f) in enumerate(mesh) - verts[NV+i] = cartesian(center(chart(mesh, f))) + verts[i] = vertices(mesh)[i] end # Now we create the faces using the new points as well # example: - # oringinal points_index (1,2,3,4) , original faces ((1,2), (2,3), (3,4),(4,5)) + # oringinal points_index (1,2,3,4) , original faces ((1,2), (2,3), (3,4),(4,1)) # new points_index(1,2,3,4,5,6,7,8), new_faces((1,5),(5,2),(2,6),(6,3),(3,7),(7,4),(4,8),(8,1)) # a,b For E=1 -> (a,c),(c,b) and so on .. - edges = similar(Edges.faces, ne) - for E in 1:NE + edges = Array{celltype(Edges)}(undef, 2NE) + for (E,p) in enumerate(Edges) + # the second half of the vertex array holds the refinement points, + # one at the center of each coarse segment + verts[NV+E] = cartesian(center(chart(Edges, p))) c = NV + E - a = Edges.faces[E][1] - b = Edges.faces[E][2] + a, b = indices(Edges, p) edges[2(E-1) + 1] = SimplexGraph(a,c) edges[2(E-1) + 2] = SimplexGraph(c,b) end + refmesh = Mesh(verts, edges) - # return the new mesh after refienments - return Mesh(verts, edges) + # return the refinement tagged with its parent so that refines recognises it + return BarycentricRefinement{U,2,T,typeof(refmesh),typeof(mesh)}(refmesh, mesh) end function barycentric_refinement(mesh::AbstractMesh{U,3}; sort=:spacefillingcurve) where U @@ -316,6 +307,8 @@ end children(mesh::Mesh{3,4}, cell) = [24*(cell-1)+1 : 24*cell] parent(mesh::BarycentricRefinement{3,4}, cell) = div((cell-1),24)+1 parent(mesh::BarycentricRefinement{U,3} where {U}, cell) = div((cell-1),6)+1 +parent(mesh::BarycentricRefinement{U,2} where {U}, cell) = div((cell-1),2)+1 + """ bisecting_refinement(mesh) -> refinement diff --git a/test/test_baryref.jl b/test/test_baryref.jl index 4645f1b..a2f410f 100644 --- a/test/test_baryref.jl +++ b/test/test_baryref.jl @@ -39,4 +39,31 @@ for T in [Float32, Float64] local b1 = skeleton(b,1) @test numcells(b1) == 2*numcells(m1) + 3*numcells(m) @test numvertices(b) == numcells(m0) + numcells(m1) -end \ No newline at end of file + + # the refinement of a curve is a BarycentricRefinement, as for surfaces + @test CompScienceMeshes.refines(Γ2, Γ) + @test numcells(Γ2) == 2*numcells(Γ) + @test numvertices(Γ2) == numvertices(Γ) + numcells(Γ) + + # parent and children are inverse of each other + local pm = CompScienceMeshes.parent(Γ2) + for E in 1:numcells(Γ) + local kids = CompScienceMeshes.children(pm, E) + @test kids == [2*(E-1)+1, 2*(E-1)+2] + @test all(CompScienceMeshes.parent(Γ2, k) == E for k in kids) + end + # curves in a three-dimensional universe, and with vertex numbering + # inherited from a parent surface + for m in (meshsegment(T(1.0), T(1)/4, 3), boundary(meshrectangle(T(1.0), T(1.0), T(1)/2, 3))) + local f = barycentric_refinement(m) + @test CompScienceMeshes.refines(f, m) + @test numcells(f) == 2*numcells(m) + @test numvertices(f) == numvertices(m) + numcells(m) + end + # a submesh of a curve refines too, which the Mesh-only signature could not do + local Λ = meshsegment(T(1.0), T(1)/4, 3) + local s = submesh((m,p) -> cartesian(CompScienceMeshes.center(chart(m,p)))[1] < T(1)/2, Λ) + local sf = barycentric_refinement(s) + @test CompScienceMeshes.refines(sf, s) + @test numcells(sf) == 2*numcells(s) +end