create grid with suitable numbering

this is to make newest vertex bisection nice-behaving

src/mesh.jl

@@ -72,7 +72,7 @@ vertices(mesh::Mesh, cell) = ntuple(i -> mesh.cells[i, cell], nvertices_cell(mes
 #    end
 #end
 
-function init_grid(m::Int, n::Int = m, v0 = (0., 0.), v1 = (1., 1.))
+function init_grid_old(m::Int, n::Int = m, v0 = (0., 0.), v1 = (1., 1.))
     r1 = LinRange(v0[1], v1[1], m + 1)
     r2 = LinRange(v0[2], v1[2], n + 1)
     coords = collect(Iterators.product(r1, r2))
@@ -92,6 +92,49 @@ function init_grid(m::Int, n::Int = m, v0 = (0., 0.), v1 = (1., 1.))
     return Mesh(vertices, cells)
 end
 
+function init_grid(m::Int, n::Int = m, v0 = (0., 0.), v1 = (1., 1.))
+    r1 = LinRange(v0[1], v1[1], m + 1)
+    r2 = LinRange(v0[2], v1[2], n + 1)
+    coords = collect(Iterators.product(r1, r2))
+
+    vertices = [x[i] for i in 1:2, x in vec(coords)]
+    cells = Array{Int, 2}(undef, 3, 2 * m * n)
+
+    vidx = zeros(Int, m + 1, n + 1)
+    k = 0
+    for j = 1 : n + 1
+	for i = 1 : m + 1
+	    if iseven(i + j)
+		vidx[i, j] = k += 1
+	    end
+	end
+    end
+    for j = 1 : n + 1
+	for i = 1 : m + 1
+	    if isodd(i + j)
+		vidx[i, j] = k += 1
+	    end
+	end
+    end
+    vidxinv = reshape(invperm(vec(vidx)), m + 1, n + 1)
+    vertices = reshape(vertices[:, vidxinv], 2, :)
+
+    e1 = CartesianIndex(1, 0)
+    e2 = CartesianIndex(0, 1)
+    k = 0
+    for I in CartesianIndices((m, n))
+	if iseven(I[1] + I[2])
+	    cells[:, k += 1] .= (vidx[I], vidx[I + e1], vidx[I + e1 + e2])
+	    cells[:, k += 1] .= (vidx[I], vidx[I + e1 + e2], vidx[I + e2])
+	else
+	    cells[:, k += 1] .= (vidx[I], vidx[I + e1], vidx[I + e2])
+	    cells[:, k += 1] .= (vidx[I + e1], vidx[I + e1 + e2], vidx[I + e2])
+	end
+    end
+
+    return Mesh(vertices, cells)
+end
+
 init_grid(img::Array{<:Any, 2}; type=:vertex) =
     type == :vertex ?
 	init_grid(size(img, 1) - 1, size(img, 2) - 1, (1.0, 1.0), size(img)) :