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/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_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 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