Optimize topological relations (#1006)

* Refactor refinement algorithms to support tuples

* Return tuples in Boundary relation

* Optimize Adjacency

* More adjustments

* Return tuple in Adjacency relation

* Optimize Coboundary

* Return tuple in Coboundary

---------

Co-authored-by: Elias Carvalho <[email protected]>

Author: juliohm

ext/mesh.jl

@@ -163,7 +163,9 @@ function vizmesh!(plot, ::Type{<:𝔼}, ::Val{2}, ::Val)
       # interleaved with a sentinel index
       inds = Int[]
       for i in 1:nfacets(t)
-        append!(inds, ∂(i))
+        for j in ∂(i)
+          push!(inds, j)
+        end
         push!(inds, nvert + 1)
       end
 

src/matrices.jl

@@ -17,8 +17,10 @@ Optionally, specify the `kind` of discretization.
 
 ## Available discretizations
 
-* `:uniform`   - `Lᵢⱼ = 1 / |𝒩(i)|, ∀j ∈ 𝒩(i)`
-* `:cotangent` - `Lᵢⱼ = cot(αᵢⱼ) + cot(βᵢⱼ), ∀j ∈ 𝒩(i)`
+* `:uniform`   - `Lᵢⱼ = 1 / |𝒜(i)|, ∀j ∈ 𝒜(i)`
+* `:cotangent` - `Lᵢⱼ = cot(αᵢⱼ) + cot(βᵢⱼ), ∀j ∈ 𝒜(i)`
+
+where `𝒜(i)` is the adjacency relation at vertex `i`.
 
 ## References
 
@@ -38,46 +40,49 @@ function laplacematrix(mesh; kind=nothing)
   𝒯 = adjusttopo(topology(mesh))
 
   # retrieve adjacency relation
-  𝒩 = Adjacency{0}(𝒯)
+  𝒜 = Adjacency{0}(𝒯)
 
   # initialize matrix
   n = nvertices(mesh)
   L = spzeros(n, n)
 
   # fill matrix
   if 𝒦 == :uniform
-    uniformlaplacian!(L, 𝒩)
+    uniformlaplacian!(L, 𝒜)
   elseif 𝒦 == :cotangent
-    cotangentlaplacian!(L, 𝒩, vertices(mesh))
+    cotangentlaplacian!(L, 𝒜, vertices(mesh))
   end
 
   L
 end
 
-function uniformlaplacian!(L, 𝒩)
+function uniformlaplacian!(L, 𝒜)
   n = size(L, 1)
   for i in 1:n
-    js = 𝒩(i)
+    js = 𝒜(i)
     for j in js
       L[i, j] = 1 / length(js)
     end
     L[i, i] = -1
   end
 end
 
-function cotangentlaplacian!(L, 𝒩, v)
+function cotangentlaplacian!(L, 𝒜, v)
   n = size(L, 1)
   for i in 1:n
-    js = CircularVector(𝒩(i))
-    for k in 1:length(js)
-      j₋, j, j₊ = js[k - 1], js[k], js[k + 1]
+    js = 𝒜(i)
+    m = length(js)
+    for k in 1:m
+      j₋ = js[mod1(k - 1, m)]
+      j = js[mod1(k, m)]
+      j₊ = js[mod1(k + 1, m)]
       vᵢ, vⱼ = v[i], v[j]
       v₋, v₊ = v[j₋], v[j₊]
       αᵢⱼ = ∠(vⱼ, v₋, vᵢ)
       βᵢⱼ = ∠(vᵢ, v₊, vⱼ)
       L[i, j] = cot(αᵢⱼ) + cot(βᵢⱼ)
     end
-    L[i, i] = -sum(L[i, js])
+    L[i, i] = -sum(j -> L[i, j], js)
   end
 end
 
@@ -100,7 +105,7 @@ function measurematrix(mesh)
   D = paramdim(mesh)
 
   # retrieve coboundary relation
-  ∂ = Coboundary{0,D}(𝒯)
+  𝒞 = Coboundary{0,D}(𝒯)
 
   # pre-compute all measures
   A = measure.(mesh)
@@ -111,7 +116,8 @@ function measurematrix(mesh)
 
   # fill matrix
   for i in 1:n
-    Aᵢ = sum(A[∂(i)]) / 3
+    js = 𝒞(i)
+    Aᵢ = sum(j -> A[j], js) / 3
     M[i, i] = 2Aᵢ
   end
 

src/refinement/catmullclark.jl

@@ -32,20 +32,20 @@ function refine(mesh, ::CatmullClark)
   # add centroids of elements
   ∂₂₀ = Boundary{2,0}(t)
   epts = map(1:nelements(t)) do elem
-    ps = view(points, ∂₂₀(elem))
-    cₒ = sum(to, ps) / length(ps)
+    is = ∂₂₀(elem)
+    cₒ = sum(j -> to(points[j]), is) / length(is)
     withcrs(mesh, cₒ)
   end
 
   # add midpoints of edges
   ∂₁₂ = Coboundary{1,2}(t)
   ∂₁₀ = Boundary{1,0}(t)
   fpts = map(1:nfacets(t)) do edge
-    ps = view(epts, ∂₁₂(edge))
-    qs = view(points, ∂₁₀(edge))
-    ∑p = sum(to, ps)
-    ∑q = sum(to, qs)
-    M = length(ps) + length(qs)
+    is = ∂₁₂(edge)
+    js = ∂₁₀(edge)
+    ∑p = sum(i -> to(epts[i]), is)
+    ∑q = sum(j -> to(points[j]), js)
+    M = length(is) + length(js)
     withcrs(mesh, (∑p + ∑q) / M)
   end
 
@@ -57,16 +57,13 @@ function refine(mesh, ::CatmullClark)
     P = to(points[u])
 
     # average of centroids
-    ps = view(epts, ∂₀₂(u))
-    F = sum(to, ps) / length(ps)
+    is = ∂₀₂(u)
+    F = sum(i -> to(epts[i]), is) / length(is)
 
     # average of midpoints
     vs = ∂₀₀(u)
     n = length(vs)
-    R = sum(vs) do v
-      uv = view(points, [u, v])
-      sum(to, uv) / 2
-    end / n
+    R = sum(v -> to(points[u]) + to(points[v]), vs) / 2n
 
     withcrs(mesh, (F + 2R + (n - 3)P) / n)
   end
@@ -81,13 +78,15 @@ function refine(mesh, ::CatmullClark)
   ∂₂₁ = Boundary{2,1}(t)
   newconnec = Connectivity{Quadrangle,4}[]
   for elem in 1:nelements(t)
-    verts = CircularVector(∂₂₀(elem))
-    edges = CircularVector(∂₂₁(elem))
-    for i in 1:length(edges)
+    verts = ∂₂₀(elem)
+    edges = ∂₂₁(elem)
+    nv = length(verts)
+    ne = length(edges)
+    for i in 1:ne
       u = elem + offset₁
-      v = edges[i] + offset₂
-      w = verts[i + 1]
-      z = edges[i + 1] + offset₂
+      v = edges[mod1(i, ne)] + offset₂
+      w = verts[mod1(i + 1, nv)]
+      z = edges[mod1(i + 1, ne)] + offset₂
       quad = connect((u, v, w, z))
       push!(newconnec, quad)
     end

src/refinement/quad.jl

@@ -21,16 +21,16 @@ function refine(mesh, ::QuadRefinement)
   # add centroids of elements
   ∂₂₀ = Boundary{2,0}(t)
   epts = map(1:nelements(t)) do elem
-    ps = view(points, ∂₂₀(elem))
-    cₒ = sum(to, ps) / length(ps)
+    is = ∂₂₀(elem)
+    cₒ = sum(j -> to(points[j]), is) / length(is)
     withcrs(mesh, cₒ)
   end
 
   # add midpoints of edges
   ∂₁₀ = Boundary{1,0}(t)
   fpts = map(1:nfacets(t)) do edge
-    ps = view(points, ∂₁₀(edge))
-    cₒ = sum(to, ps) / length(ps)
+    is = ∂₁₀(edge)
+    cₒ = sum(i -> to(points[i]), is) / length(is)
     withcrs(mesh, cₒ)
   end
 
@@ -47,13 +47,15 @@ function refine(mesh, ::QuadRefinement)
   ∂₂₁ = Boundary{2,1}(t)
   newconnec = Connectivity{Quadrangle,4}[]
   for elem in 1:nelements(t)
-    verts = CircularVector(∂₂₀(elem))
-    edges = CircularVector(∂₂₁(elem))
-    for i in 1:length(edges)
+    verts = ∂₂₀(elem)
+    edges = ∂₂₁(elem)
+    nv = length(verts)
+    ne = length(edges)
+    for i in 1:ne
       u = elem + offset₁
-      v = edges[i] + offset₂
-      w = verts[i + 1]
-      z = edges[i + 1] + offset₂
+      v = edges[mod1(i, ne)] + offset₂
+      w = verts[mod1(i + 1, nv)]
+      z = edges[mod1(i + 1, ne)] + offset₂
       quad = connect((u, v, w, z))
       push!(newconnec, quad)
     end

src/refinement/tri.jl

@@ -24,8 +24,8 @@ function refine(mesh, ::TriRefinement)
   # add centroids of elements
   ∂₂₀ = Boundary{2,0}(t)
   epts = map(1:nelements(t)) do elem
-    ps = view(points, ∂₂₀(elem))
-    cₒ = sum(to, ps) / length(ps)
+    is = ∂₂₀(elem)
+    cₒ = sum(i -> to(points[i]), is) / length(is)
     withcrs(mesh, cₒ)
   end
 
@@ -40,11 +40,12 @@ function refine(mesh, ::TriRefinement)
   # connect vertices into new triangles
   newconnec = Connectivity{Triangle,3}[]
   for elem in 1:nelements(t)
-    verts = CircularVector(∂₂₀(elem))
-    for i in 1:length(verts)
+    verts = ∂₂₀(elem)
+    nv = length(verts)
+    for i in 1:nv
       u = elem + offset
-      v = verts[i]
-      w = verts[i + 1]
+      v = verts[mod1(i, nv)]
+      w = verts[mod1(i + 1, nv)]
       tri = connect((u, v, w))
       push!(newconnec, tri)
     end

src/refinement/trisubdivision.jl

@@ -36,10 +36,10 @@ function refine(mesh, ::TriSubdivision)
   midpoints = Dict{Tuple{Int,Int},Int}()
   ∂₁₀ = Boundary{1,0}(t)
   for eind in 1:nfacets(t)
-    i, j = sort(∂₁₀(eind))
+    i, j = ∂₁₀(eind)
     edge = Segment(points[i], points[j])
     push!(points, centroid(edge))
-    midpoints[(i, j)] = (np += 1)
+    midpoints[_ordered(i, j)] = (np += 1)
   end
 
   # construct subtriangles of faces

src/topologies/grid.jl

@@ -120,15 +120,15 @@ nvertices(t::GridTopology) = prod(t.dims .+ t.open)
 function element(t::GridTopology{D}, ind) where {D}
   ∂ = Boundary{D,0}(t)
   T = elementtype(t)
-  connect(Tuple(∂(ind)), T)
+  connect(∂(ind), T)
 end
 
 nelements(t::GridTopology) = prod(t.dims)
 
 function facet(t::GridTopology{D}, ind) where {D}
   ∂ = Boundary{D - 1,0}(t)
   T = facettype(t)
-  connect(Tuple(∂(ind)), T)
+  connect(∂(ind), T)
 end
 
 nfacets(t::GridTopology{1}) = t.dims[1] + t.open[1]

src/toporelations/adjacency.jl

@@ -34,6 +34,7 @@ function (𝒜::Adjacency{0,D,T})(ind::Integer) where {D,T<:GridTopology}
   # construct topology for vertices
   vtopo = GridTopology(dims .+ 1, cycl)
   𝒜vert = Adjacency{D}(vtopo)
+
   𝒜vert(ind)
 end
 
@@ -44,25 +45,25 @@ function (𝒜::Adjacency{D,D,T})(ind::Integer) where {D,T<:GridTopology}
   cycl = isperiodic(topo)
   cind = elem2cart(topo, ind)
 
-  # offsets along each dimension
-  offsets = [ntuple(i -> i == d ? s : 0, D) for d in 1:D for s in (-1, 1)]
+  inds = Int[]
+  for d in 1:D, s in (-1, 1)
+    # offset along each dimension
+    offset = ntuple(i -> i == d ? s : 0, D)
 
-  ninds = NTuple{D,Int}[]
-  for offset in offsets
     # apply offset to center index
     sind = cind .+ offset
 
     # wrap indices in case of periodic dimension
-    wrap(i) = mod1(sind[i], dims[i])
-    wind = ntuple(i -> cycl[i] ? wrap(i) : sind[i], D)
+    wind = ntuple(D) do i
+      cycl[i] ? mod1(sind[i], dims[i]) : sind[i]
+    end
 
     # discard invalid indices
     valid(i) = 1 ≤ wind[i] ≤ dims[i]
-    all(valid, 1:D) && push!(ninds, wind)
+    all(valid, 1:D) && push!(inds, cart2elem(topo, wind...))
   end
 
-  # return linear index of element
-  [cart2elem(topo, ind...) for ind in ninds]
+  ntuple(i -> inds[i], length(inds))
 end
 
 # -------------------
@@ -74,32 +75,32 @@ function (𝒜::Adjacency{0,2,T})(vert::Integer) where {T<:HalfEdgeTopology}
   e = half4vert(𝒜.topology, vert)
 
   # initialize result
-  vertices = [e.half.head]
+  inds = [e.half.head]
 
   # search in CCW orientation
   p = e.prev
   h = p.half
   while !isnothing(h.elem) && h != e
-    push!(vertices, p.head)
+    push!(inds, p.head)
     p = h.prev
     h = p.half
   end
 
   # if border edge is hit
   if isnothing(h.elem)
     # add last arm manually
-    push!(vertices, p.head)
+    push!(inds, p.head)
 
     # search in CW orientation
     h = e.half
     while !isnothing(h.elem)
       n = h.next
       h = n.half
-      pushfirst!(vertices, h.head)
+      pushfirst!(inds, h.head)
     end
   end
 
-  vertices
+  ntuple(i -> inds[i], length(inds))
 end
 
 # adjacent elements in a 2D half-edge topology
@@ -117,5 +118,5 @@ function (𝒜::Adjacency{2,2,T})(ind::Integer) where {T<:HalfEdgeTopology}
     n = n.next
   end
 
-  inds
+  ntuple(i -> inds[i], length(inds))
 end