Fix adjacency matrix computation if mesh contains undead vertices.

src/Triangle.Tests/Tools/AdjacencyMatrixTest.cs

@@ -40,20 +40,19 @@ namespace TriangleNet.Tests.Tools
 
             var ai = matrix.RowIndices;
 
-            // Highest vertex id after renumbering is 5, since the duplicate
-            // vertex is still present in the mesh vertices list.
-            Assert.AreEqual(5, ai.Max());
+            // Highest vertex id after renumbering is 4, duplicates
+            // are ignored.
+            Assert.AreEqual(4, ai.Max());
 
             // Get the single, duplicate vertex.
             var dup = mesh.Vertices
                 .Where(v => v.Type == VertexType.UndeadVertex)
                 .Single();
 
-            // The duplicate vertex is part of the matrix, since it is assumed
-            // to be adjacent to itself.
-            Assert.AreEqual(1, ai.Count(i => i == dup.id));
-
-            // TODO: fix AdjacencyMatrix!!!
+            // Side effect: undead vertices will have negative indices
+            // after computing the adjacency matrix.
+            Assert.IsTrue(dup.id < 0);
+            Assert.IsTrue(!ai.Contains(dup.id));
         }
 
         private List<Vertex> GetVertices(bool includeDuplicate)

src/Triangle/Tools/AdjacencyMatrix.cs

@@ -28,38 +28,46 @@ namespace TriangleNet.Tools
         /// <summary>
         /// Gets the number of columns (nodes of the mesh).
         /// </summary>
-        public readonly int N;
+        public readonly int ColumnCount;
 
         /// <summary>
         /// Gets the column pointers.
         /// </summary>
-        public int[] ColumnPointers
-        {
-            get { return pcol; }
-        }
+        public int[] ColumnPointers => pcol;
 
         /// <summary>
         /// Gets the row indices.
         /// </summary>
-        public int[] RowIndices
-        {
-            get { return irow; }
-        }
+        public int[] RowIndices => irow;
 
         /// <summary>
         /// Initializes a new instance of the <see cref="AdjacencyMatrix" /> class.
         /// </summary>
         /// <param name="mesh">The mesh.</param>
+        /// <remarks>
+        /// NOTE: as a side effect, computing the adjacency matrix will affect the
+        /// node numbering of the mesh.
+        /// </remarks>
         public AdjacencyMatrix(Mesh mesh)
         {
-            this.N = mesh.vertices.Count;
+            int n = mesh.vertices.Count;
+
+            // Undead vertices should not be considered in the adjacency matrix.
+            ColumnCount = n - mesh.undeads;
+
+            // Renumber nodes, excluding undeads.
+            int i = 0;
+            foreach (var vertex in mesh.vertices.Values)
+            {
+                vertex.id = vertex.type == VertexType.UndeadVertex ? -i : i++;
+            }
 
             // Set up the adj_row adjacency pointer array.
-            this.pcol = AdjacencyCount(mesh);
-            this.nnz = pcol[N];
+            pcol = AdjacencyCount(mesh);
+            nnz = pcol[ColumnCount];
 
             // Set up the adj adjacency array.
-            this.irow = AdjacencySet(mesh, this.pcol);
+            irow = AdjacencySet(mesh, pcol);
 
             SortIndices();
         }
@@ -72,16 +80,16 @@ namespace TriangleNet.Tools
         /// <exception cref="ArgumentException"></exception>
         public AdjacencyMatrix(int[] pcol, int[] irow)
         {
-            this.N = pcol.Length - 1;
+            ColumnCount = pcol.Length - 1;
 
-            this.nnz = pcol[N];
+            nnz = pcol[ColumnCount];
 
             this.pcol = pcol;
             this.irow = irow;
 
             if (pcol[0] != 0)
             {
-                throw new ArgumentException("Expected 0-based indexing.", "pcol");
+                throw new ArgumentException("Expected 0-based indexing.", nameof(pcol));
             }
 
             if (irow.Length < nnz)
@@ -104,7 +112,7 @@ namespace TriangleNet.Tools
             band_lo = 0;
             band_hi = 0;
 
-            for (i = 0; i < N; i++)
+            for (i = 0; i < ColumnCount; i++)
             {
                 for (j = pcol[i]; j < pcol[i + 1]; j++)
                 {
@@ -134,7 +142,7 @@ namespace TriangleNet.Tools
         /// </remarks>
         int[] AdjacencyCount(Mesh mesh)
         {
-            int n = N;
+            int n = ColumnCount;
             int n1, n2, n3;
             int tid, nid;
 
@@ -211,7 +219,7 @@ namespace TriangleNet.Tools
         /// </remarks>
         int[] AdjacencySet(Mesh mesh, int[] pcol)
         {
-            int n = this.N;
+            int n = ColumnCount;
 
             int[] col = new int[n];
 
@@ -280,9 +288,9 @@ namespace TriangleNet.Tools
         /// </summary>
         public void SortIndices()
         {
-            int k1, k2, n = N;
+            int k1, k2, n = ColumnCount;
 
-            int[] list = this.irow;
+            var list = irow;
 
             // Ascending sort the entries for each column.
             for (int i = 0; i < n; i++)

src/Triangle/Tools/CuthillMcKee.cs

@@ -25,9 +25,6 @@ namespace TriangleNet.Tools
         /// <returns>Permutation vector.</returns>
         public int[] Renumber(Mesh mesh)
         {
-            // Algorithm needs linear numbering of the nodes.
-            mesh.Renumber(NodeNumbering.Linear);
-
             return Renumber(new AdjacencyMatrix(mesh));
         }
 
@@ -58,7 +55,7 @@ namespace TriangleNet.Tools
 
             if (Log.Verbose)
             {
-                Log.Instance.Info(String.Format("Reverse Cuthill-McKee (Bandwidth: {0} > {1})",
+                Log.Instance.Info(string.Format("Reverse Cuthill-McKee (Bandwidth: {0} > {1})",
                     bandwidth1, bandwidth2));
             }
 
@@ -81,7 +78,7 @@ namespace TriangleNet.Tools
         int[] GenerateRcm()
         {
             // Number of nodes in the mesh.
-            int n = matrix.N;
+            int n = matrix.ColumnCount;
 
             int[] perm = new int[n];
 
@@ -167,7 +164,7 @@ namespace TriangleNet.Tools
             int nbr, node;
 
             // Number of nodes in the mesh.
-            int n = matrix.N;
+            int n = matrix.ColumnCount;
 
             // Workspace, int DEG[NODE_NUM], a temporary vector used to hold 
             // the degree of the nodes in the section graph specified by mask and root.
@@ -603,7 +600,7 @@ namespace TriangleNet.Tools
             int band_lo = 0;
             int band_hi = 0;
 
-            int n = matrix.N;
+            int n = matrix.ColumnCount;
 
             for (i = 0; i < n; i++)
             {
@@ -625,7 +622,7 @@ namespace TriangleNet.Tools
         /// <returns>The inverse permutation.</returns>
         int[] PermInverse(int[] perm)
         {
-            int n = matrix.N;
+            int n = matrix.ColumnCount;
 
             int[] perm_inv = new int[n];