Refactor adjsort into adjsortperm

src/topologies/halfedge.jl

@@ -144,9 +144,8 @@
 
   # sort elements to make sure that they
   # are traversed in adjacent-first order
-  elemsort = sort ? adjsort(elems) : elems
-  eleminds = sort ? indexin(elemsort, elems) : 1:length(elems)
-  adjelems::Vector{Vector{Int}} = map(collect ∘ indices, elemsort)
+  eleminds = sort ? adjsortperm(elems) : eachindex(elems)
+  adjelems::Vector{Vector{Int}} = map(collect ∘ indices ∘ Base.Fix1(getindex, elems), eleminds)
 
   # start assuming that all elements are
   # oriented consistently (e.g. CCW)
@@ -178,7 +177,7 @@
       elem = eleminds[other]
       inds = adjelems[other]
       n = length(inds)
-      if anyhalf(inds, half4pair) || disconnected
+      if anyhalf(half4pair, inds) || disconnected
         # at least one half-edge has been found, so we can assess
         # the orientation w.r.t. previously added elements/edges
         added = true
@@ -191,7 +190,7 @@
 
       # if the other element is already claimed by any half-edge
       # then the element must be reversed before further updates
-      isreversed[other] = anyhalfclaimed(inds, half4pair)
+      isreversed[other] = anyhalfclaimed(half4pair, inds)
 
       # insert half-edges in consistent orientation
       if isreversed[other]
@@ -332,45 +331,87 @@
 # HELPER FUNCTIONS
 # -----------------
 
-function adjsort(elems::AbstractVector{<:Connectivity})
-  # initialize list of adjacent elements
-  # with first element from original list
-  list = indices.(elems)
-  adjs = Tuple[popfirst!(list)]
-
-  # the loop will terminate if the mesh
-  # is manifold, and that is always true
-  # with half-edge topology
-  while !isempty(list)
-    # lookup all elements that share at least
-    # one vertex with the last adjacent element
-    found = false
-    vinds = last(adjs)
-    for i in vinds
-      einds = findall(e -> i ∈ e, list)
-      if !isempty(einds)
-        # lookup all elements that share at
-        # least two vertices (i.e. edge)
-        for j in sort(einds, rev=true)
-          if length(vinds ∩ list[j]) > 1
-            found = true
-            push!(adjs, popat!(list, j))
-          end
-        end
+# permutation of elements in adjacent-first order
+function adjsortperm(elems::AbstractVector{<:Connectivity})
+  reduce(vcat, conneccomps(elems))
+end
+
+# connected components from list of elements
+function conneccomps(elems::AbstractVector{<:Connectivity})
+  # initialize list of connected components
+  comps = [[firstindex(elems)]]
+
+  # initialize list of seen vertices
+  seen = Set{Int}()
+  for v in indices(first(elems))
+    push!(seen, v)
+  end
+
+  # remaining elements to process
+  remaining = collect(eachindex(elems)[2:end])
+
+  added = false
+  while !isempty(remaining)
+    iter = 1
+    while iter ≤ length(remaining)
+      elem = elems[remaining[iter]]
+
+      # manually union-split two most common polytopes
+      # for type stability and maximum performance
+      isadjacent = if elem isa Connectivity{Triangle,3}
+        adjelem!(seen, elem)
+      elseif elem isa Connectivity{Quadrangle,4}
+        adjelem!(seen, elem)
+      else
+        adjelem!(seen, elem)
+      end
+
+      if isadjacent
+        push!(last(comps), popat!(remaining, iter))
+        added = true
+      else
+        iter += 1
       end
     end
 
-    if !found && !isempty(list)
-      # we are done with this connected component
-      # pop a new element from the original list
-      push!(adjs, popfirst!(list))
+    if added
+      # new vertices were "seen" while iterating `remaining`, so
+      # we need to iterate again because there may be elements
+      # which are now adjacent with the newly "seen" vertices
+      added = false
+    elseif !isempty(remaining)
+      # there are more elements, but none are adjacent to
+      # previously seen elements; pop a new element from
+      # the original list to start a new connected component
+      push!(comps, Int[])
+      push!(last(comps), popfirst!(remaining))
+
+      # a disconnected component means that ≥n-1 vertices in
+      # the newest element haven't been "seen"; its possible
+      # the new component is connected by a single vertex
+      for v in indices(elems[last(last(comps))])
+        push!(seen, v)
+      end
     end
   end
 
-  connect.(adjs)
+  comps
+end
+
+# update seen vertices if there are ≥2 common indices
+function adjelem!(seen, elem)
+  inds = indices(elem)
+  if count(∈(seen), inds) > 1
+    for v in inds
+      push!(seen, v)
+    end
+    return true
+  end
+  return false
 end
 
-function anyhalf(inds, half4pair)
+# true if the half-edges already contain the indices
+function anyhalf(half4pair, inds)
   n = length(inds)
   for i in eachindex(inds)
     uv = (inds[i], inds[mod1(i + 1, n)])
@@ -381,7 +422,8 @@
   return false
 end
 
-function anyhalfclaimed(inds, half4pair)
+# true if the half-edges already have elements assigned to the indices
+function anyhalfclaimed(half4pair, inds)
   n = length(inds)
   for i in eachindex(inds)
     uv = (inds[i], inds[mod1(i + 1, n)])
@@ -392,5 +434,6 @@
   return false
 end
 
+# integer addition mod1
 add0(i, n) = i
 add1(i, n) = mod1(i + 1, n)

test/topologies.jl

@@ -512,8 +512,8 @@ end
   g = GridTopology(10, 10)
   t = convert(HalfEdgeTopology, g)
   @test t[begin] == connect((1, 2, 13, 12), Quadrangle)
-  @test t[end] == connect((120, 109, 110, 121), Quadrangle)
-  @test t[10] == connect((22, 21, 10, 11), Quadrangle)
+  @test t[end] == connect((109, 110, 121, 120), Quadrangle)
+  @test t[10] == connect((21, 10, 11, 22), Quadrangle)
   @test length(t) == 100
   @test eltype(t) == Connectivity{Quadrangle,4}
   for e in t

test/toporelations.jl

@@ -488,26 +488,26 @@ end
   @test ∂(4) == (1, 5, 3)
   ∂ = Boundary{2,1}(t)
   @test ∂(1) == (1, 2, 3, 4)
-  @test ∂(2) == (2, 7, 8)
-  @test ∂(3) == (3, 8, 9, 5)
-  @test ∂(4) == (4, 5, 6)
+  @test ∂(2) == (2, 5, 6)
+  @test ∂(3) == (3, 6, 7, 8)
+  @test ∂(4) == (4, 8, 9)
   ∂ = Boundary{1,0}(t)
   @test ∂(1) == (1, 2)
   @test ∂(2) == (2, 6)
   @test ∂(3) == (6, 5)
   @test ∂(4) == (5, 1)
-  @test ∂(5) == (5, 3)
-  @test ∂(6) == (3, 1)
-  @test ∂(7) == (2, 4)
-  @test ∂(8) == (4, 6)
-  @test ∂(9) == (4, 3)
+  @test ∂(5) == (2, 4)
+  @test ∂(6) == (4, 6)
+  @test ∂(7) == (4, 3)
+  @test ∂(8) == (3, 5)
+  @test ∂(9) == (3, 1)
   𝒞 = Coboundary{0,1}(t)
-  @test 𝒞(1) == (1, 4, 6)
-  @test 𝒞(2) == (7, 2, 1)
-  @test 𝒞(3) == (6, 5, 9)
-  @test 𝒞(4) == (9, 8, 7)
-  @test 𝒞(5) == (3, 5, 4)
-  @test 𝒞(6) == (2, 8, 3)
+  @test 𝒞(1) == (1, 4, 9)
+  @test 𝒞(2) == (5, 2, 1)
+  @test 𝒞(3) == (9, 8, 7)
+  @test 𝒞(4) == (7, 6, 5)
+  @test 𝒞(5) == (3, 8, 4)
+  @test 𝒞(6) == (2, 6, 3)
   𝒞 = Coboundary{0,2}(t)
   @test 𝒞(1) == (1, 4)
   @test 𝒞(2) == (2, 1)
@@ -520,11 +520,11 @@ end
   @test 𝒞(2) == (1, 2)
   @test 𝒞(3) == (1, 3)
   @test 𝒞(4) == (1, 4)
-  @test 𝒞(5) == (4, 3)
-  @test 𝒞(6) == (4,)
-  @test 𝒞(7) == (2,)
-  @test 𝒞(8) == (2, 3)
-  @test 𝒞(9) == (3,)
+  @test 𝒞(5) == (2,)
+  @test 𝒞(6) == (2, 3)
+  @test 𝒞(7) == (3,)
+  @test 𝒞(8) == (3, 4)
+  @test 𝒞(9) == (4,)
   𝒜 = Adjacency{0}(t)
   @test 𝒜(1) == (2, 5, 3)
   @test 𝒜(2) == (4, 6, 1)