Dof reindexing for FESpaces

NEWS.md

@@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+
+- Added `reindex_free_dof_ids(space, algorithm)` to reorder the free DOFs of a `FESpace` for bandwidth/profile/fill-in reduction. Supported algorithms: `:rcm` (Reverse Cuthill-McKee), `:sloan` (wavefront minimisation), and `:coordinates` (sort by spatial coordinates, for Lagrangian spaces). An externally computed permutation vector can also be applied directly via `reindex_free_dof_ids(space, free_dof_ids)`. Since PR[#1299](https://github.com/gridap/Gridap.jl/pull/1299).
+
 ### Changed
 
 - Minor generalisation of `return_value` for `LinearCombinationMap`. Since PR[#1298](https://github.com/gridap/Gridap.jl/pull/1298).

Project.toml

@@ -1,7 +1,7 @@
 name = "Gridap"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
-authors = ["Santiago Badia <santiago.badia@monash.edu>", "Francesc Verdugo <f.verdugo.rojano@vu.nl>", "Alberto F. Martin <alberto.f.martin@anu.edu.au>"]
 version = "0.20.6"
+authors = ["Santiago Badia <santiago.badia@monash.edu>", "Francesc Verdugo <f.verdugo.rojano@vu.nl>", "Alberto F. Martin <alberto.f.martin@anu.edu.au>"]
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -40,6 +40,7 @@ TikzPictures = "37f6aa50-8035-52d0-81c2-5a1d08754b2d"
 TikzPicturesExt = "TikzPictures"
 
 [compat]
+AMD = "0.5"
 AbstractTrees = "0.3.3, 0.4"
 Aqua = "0.8"
 AutoHashEquals = "2.2.0"
@@ -73,8 +74,9 @@ WriteVTK = "1.21.1"
 julia = "1.10"
 
 [extras]
+AMD = "14f7f29c-3bd6-536c-9a0b-7339e30b5a3e"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Test"]
+test = ["AMD", "Aqua", "Test"]

src/Arrays/Arrays.jl

@@ -111,6 +111,7 @@ export scatter_table_values
 export gather_posneg_table_values
 export gather_posneg_table_values!
 export scatter_posneg_table_values
+export inverse_table
 
 export IdentityVector
 

src/Arrays/Tables.jl

@@ -419,7 +419,9 @@ function inverse_table(
   o = one(P)
   ptrs = zeros(P,nb+1)
   @inbounds for b in a_to_lb_to_b_data
-    ptrs[b+1] += o
+    if b > 0 
+      ptrs[b+1] += o
+    end
   end
   length_to_ptrs!(ptrs)
 
@@ -430,7 +432,7 @@ function inverse_table(
     e = a_to_lb_to_b_ptrs[a+1] - o
     @inbounds for p in s:e
       b = a_to_lb_to_b_data[p]
-      if b != UNSET
+      if b > 0
         data[ptrs[b]] = a
         ptrs[b] += o
       end
@@ -1034,3 +1036,71 @@ end
 function Base.copy(a::Table)
   Table(copy(a.data),copy(a.ptrs))
 end
+
+"""
+    compute_adjacency(t::Table[, nfree]) -> Table
+
+Build the symmetric node adjacency graph from a hyper-edge table `t`
+(e.g. a cell-to-DOF table).  Two nodes are adjacent iff they appear
+in the same row.  Entries ≤ 0 are treated as inactive (e.g. Dirichlet
+DOF IDs encoded as negative values) and skipped.
+
+`nfree` defaults to the maximum positive entry in `t`.
+
+Returns `adj` where `adj[i]` lists the neighbours of node `i` in sorted
+order.  Node degrees are available cheaply as `diff(adj.ptrs)`.
+"""
+function compute_adjacency(
+  t::Table{T}, nfree::Int = Int(maximum(t.data; init=zero(T)))
+) where T
+  d2c = inverse_table(t, nfree)
+
+  # Upper-bound pointer array: degree(i) ≤ Σ |row(c)| for cells c containing i
+  ub_ptrs = Vector{Int}(undef, nfree + 1)
+  ub_ptrs[1] = 1
+  for i in 1:nfree
+    s = 0
+    for k in datarange(d2c, i)
+      s += length(datarange(t, d2c.data[k]))
+    end
+    ub_ptrs[i+1] = ub_ptrs[i] + s
+  end
+
+  # Collect positive, non-self neighbours (may contain duplicates)
+  ub_data  = Vector{Int32}(undef, ub_ptrs[nfree+1] - 1)
+  fill_end = copy(ub_ptrs)
+  for i in 1:nfree
+    for k in datarange(d2c, i)
+      c = d2c.data[k]
+      for k2 in datarange(t, c)
+        j = t.data[k2]
+        if j > 0 && j != i
+          ub_data[fill_end[i]] = j
+          fill_end[i] += 1
+        end
+      end
+    end
+  end
+
+  # Sort each row slice, then count and copy unique neighbours
+  ptrs = Vector{Int}(undef, nfree + 1)
+  ptrs[1] = 1
+  for i in 1:nfree
+    sort!(view(ub_data, ub_ptrs[i]:fill_end[i]-1))
+    n = 0; prev = zero(Int32)
+    for k in ub_ptrs[i]:fill_end[i]-1
+      v = ub_data[k]; v != prev && (n += 1; prev = v)
+    end
+    ptrs[i+1] = ptrs[i] + n
+  end
+
+  data = Vector{Int32}(undef, ptrs[nfree+1] - 1)
+  for i in 1:nfree
+    p = ptrs[i]; prev = zero(Int32)
+    for k in ub_ptrs[i]:fill_end[i]-1
+      v = ub_data[k]; v != prev && (data[p] = v; p += 1; prev = v)
+    end
+  end
+
+  return Table(data, ptrs)
+end