Bugfix: Curve dual support

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

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)
 

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
+
+    # 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

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