diff --git a/Triangle.NET/MeshRenderer.Core/RenderData.cs b/Triangle.NET/MeshRenderer.Core/RenderData.cs index c9a7084..c97dc96 100644 --- a/Triangle.NET/MeshRenderer.Core/RenderData.cs +++ b/Triangle.NET/MeshRenderer.Core/RenderData.cs @@ -102,7 +102,7 @@ namespace MeshRenderer.Core // Copy segments n = mesh.Segments.Count; - if (n > 0) + if (n > 0 && mesh.IsPolygon) { var segments = new List(2 * n); foreach (var seg in mesh.Segments) diff --git a/Triangle.NET/TestApp/FormMain.cs b/Triangle.NET/TestApp/FormMain.cs index 240964c..c51fb46 100644 --- a/Triangle.NET/TestApp/FormMain.cs +++ b/Triangle.NET/TestApp/FormMain.cs @@ -111,7 +111,7 @@ namespace MeshExplorer if (input != null) { - settings.CurrentFile = "GEN"; + settings.CurrentFile = "tmp-" + DateTime.Now.ToString("HH-mm-ss"); HandleNewInput(); } } @@ -215,6 +215,20 @@ namespace MeshExplorer #region State changes + private void LockOnException() + { + btnMesh.Enabled = false; + btnSmooth.Enabled = false; + + //menuFileSave.Enabled = false; + //menuFileExport.Enabled = false; + menuViewVoronoi.Enabled = false; + menuToolsCheck.Enabled = false; + menuToolsRcm.Enabled = false; + + settings.ExceptionThrown = true; + } + private void HandleNewInput() { // Reset mesh @@ -405,7 +419,7 @@ namespace MeshExplorer if (meshControlView.ParamQualityChecked) { btnMesh.Text = "Refine"; - btnSmooth.Enabled = true; + btnSmooth.Enabled = mesh.IsPolygon; } } else if (meshControlView.ParamQualityChecked) @@ -449,7 +463,7 @@ namespace MeshExplorer mesh.SetOption(Options.Convex, true); } - //try + try { //sw.Start(); mesh.Triangulate(input); @@ -464,11 +478,11 @@ namespace MeshExplorer settings.RefineMode = true; } } - //catch (Exception ex) - //{ - // settings.ExceptionThrown = true; - // DarkMessageBox.Show("Exception - Triangulate", ex.Message); - //} + catch (Exception ex) + { + LockOnException(); + DarkMessageBox.Show("Exception - Triangulate", ex.Message, MessageBoxButtons.OK); + } UpdateLog(); } @@ -500,8 +514,8 @@ namespace MeshExplorer } catch (Exception ex) { - settings.ExceptionThrown = true; - DarkMessageBox.Show("Exception - Refine", ex.Message); + LockOnException(); + DarkMessageBox.Show("Exception - Refine", ex.Message, MessageBoxButtons.OK); } UpdateLog(); @@ -524,6 +538,11 @@ namespace MeshExplorer { if (mesh == null || settings.ExceptionThrown) return; + if (!mesh.IsPolygon) + { + return; + } + Stopwatch sw = new Stopwatch(); try @@ -538,8 +557,8 @@ namespace MeshExplorer } catch (Exception ex) { - settings.ExceptionThrown = true; - DarkMessageBox.Show("Exception - Smooth", ex.Message); + LockOnException(); + DarkMessageBox.Show("Exception - Smooth", ex.Message, MessageBoxButtons.OK); } UpdateLog(); diff --git a/Triangle.NET/TestApp/Generators/StarInBox.cs b/Triangle.NET/TestApp/Generators/StarInBox.cs index ddfd762..ecc40ea 100644 --- a/Triangle.NET/TestApp/Generators/StarInBox.cs +++ b/Triangle.NET/TestApp/Generators/StarInBox.cs @@ -71,20 +71,20 @@ namespace MeshExplorer.Generators x = r * Math.Cos(i * step + e); y = r * Math.Sin(i * step + e); - input.AddPoint(x, y); - input.AddSegment(0, i + 1); + input.AddPoint(x, y, 2); + input.AddSegment(0, i + 1, 2); } - input.AddPoint(-1, -1); // Box - input.AddPoint(1, -1); - input.AddPoint(1, 1); - input.AddPoint(-1, 1); + input.AddPoint(-1, -1, 1); // Box + input.AddPoint(1, -1, 1); + input.AddPoint(1, 1, 1); + input.AddPoint(-1, 1, 1); numRays = input.Count; - input.AddSegment(numRays - 1, numRays - 2); - input.AddSegment(numRays - 2, numRays - 3); - input.AddSegment(numRays - 3, numRays - 4); - input.AddSegment(numRays - 4, numRays - 1); + input.AddSegment(numRays - 1, numRays - 2, 1); + input.AddSegment(numRays - 2, numRays - 3, 1); + input.AddSegment(numRays - 3, numRays - 4, 1); + input.AddSegment(numRays - 4, numRays - 1, 1); return input; } diff --git a/Triangle.NET/TestApp/Views/AboutView.Designer.cs b/Triangle.NET/TestApp/Views/AboutView.Designer.cs index 05516dc..b892f81 100644 --- a/Triangle.NET/TestApp/Views/AboutView.Designer.cs +++ b/Triangle.NET/TestApp/Views/AboutView.Designer.cs @@ -79,7 +79,7 @@ this.label19.Name = "label19"; this.label19.Size = new System.Drawing.Size(134, 40); this.label19.TabIndex = 6; - this.label19.Text = "Beta 2 (2012-10-14)\r\nChristian Woltering\r\nMIT"; + this.label19.Text = "Beta 3 (2012-10-30)\r\nChristian Woltering\r\nMIT"; // // label18 // diff --git a/Triangle.NET/Triangle/Data/Otri.cs b/Triangle.NET/Triangle/Data/Otri.cs index 06eef19..b42e252 100644 --- a/Triangle.NET/Triangle/Data/Otri.cs +++ b/Triangle.NET/Triangle/Data/Otri.cs @@ -16,7 +16,7 @@ namespace TriangleNet.Data /// An oriented triangle. /// /// - /// Iincludes a pointer to a triangle and orientation. + /// Includes a pointer to a triangle and orientation. /// The orientation denotes an edge of the triangle. Hence, there are /// three possible orientations. By convention, each edge always points /// counterclockwise about the corresponding triangle. diff --git a/Triangle.NET/Triangle/Data/Segment.cs b/Triangle.NET/Triangle/Data/Segment.cs index 183cd39..9e7263a 100644 --- a/Triangle.NET/Triangle/Data/Segment.cs +++ b/Triangle.NET/Triangle/Data/Segment.cs @@ -69,14 +69,6 @@ namespace TriangleNet.Data get { return this.vertices[1].id; } } - /// - /// Gets the segments endpoint. - /// - public Vertex GetVertex(int index) - { - return this.vertices[index]; // TODO: Check range? - } - /// /// Gets the segment boundary mark. /// @@ -85,6 +77,14 @@ namespace TriangleNet.Data get { return this.boundary; } } + /// + /// Gets the segments endpoint. + /// + public Vertex GetVertex(int index) + { + return this.vertices[index]; // TODO: Check range? + } + #endregion public override int GetHashCode() diff --git a/Triangle.NET/Triangle/Smoothing/SimpleSmoother.cs b/Triangle.NET/Triangle/Smoothing/SimpleSmoother.cs index 77e4bdd..37c7138 100644 --- a/Triangle.NET/Triangle/Smoothing/SimpleSmoother.cs +++ b/Triangle.NET/Triangle/Smoothing/SimpleSmoother.cs @@ -16,6 +16,10 @@ namespace TriangleNet.Smoothing /// /// Simple mesh smoother implementation. /// + /// + /// Vertices wich should not move (e.g. segment vertices) MUST have a + /// boundary mark greater than 0. + /// public class SimpleSmoother : ISmoother { Mesh mesh; diff --git a/Triangle.NET/Triangle/Tools/AdjacencyMatrix.cs b/Triangle.NET/Triangle/Tools/AdjacencyMatrix.cs new file mode 100644 index 0000000..f773033 --- /dev/null +++ b/Triangle.NET/Triangle/Tools/AdjacencyMatrix.cs @@ -0,0 +1,404 @@ +// ----------------------------------------------------------------------- +// +// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/ +// +// ----------------------------------------------------------------------- + +namespace TriangleNet.Tools +{ + using System; + using System.Collections.Generic; + using System.Linq; + using System.Text; + + /// + /// The adjacency matrix of the mesh. + /// + public class AdjacencyMatrix + { + // Number of nodes in the mesh. + int node_num; + + // Number of adjacency entries. + int adj_num; + + // Pointers into the actual adjacency structure adj. Information about row k is + // stored in entries adj_row(k) through adj_row(k+1)-1 of adj. Size: node_num + 1 + int[] adj_row; + + // The adjacency structure. For each row, it contains the column indices + // of the nonzero entries. Size: adj_num + int[] adj; + + public int[] AdjacencyRow + { + get { return adj_row; } + } + + public int[] Adjacency + { + get { return adj; } + } + + public AdjacencyMatrix(Mesh mesh) + { + this.node_num = mesh.vertices.Count; + + // Set up the adj_row adjacency pointer array. + this.adj_row = AdjacencyCount(mesh); + this.adj_num = adj_row[node_num] - 1; + + // Set up the adj adjacency array. + this.adj = AdjacencySet(mesh, this.adj_row); + } + + /// + /// Computes the bandwidth of an adjacency matrix. + /// + /// Bandwidth of the adjacency matrix. + public int Bandwidth() + { + int band_hi; + int band_lo; + int col; + int i, j; + + band_lo = 0; + band_hi = 0; + + for (i = 0; i < node_num; i++) + { + for (j = adj_row[i]; j <= adj_row[i + 1] - 1; j++) + { + col = adj[j - 1]; + band_lo = Math.Max(band_lo, i - col); + band_hi = Math.Max(band_hi, col - i); + } + } + + return band_lo + 1 + band_hi; + } + + #region Adjacency matrix + + /// + /// Counts adjacencies in a triangulation. + /// + /// + /// This routine is called to count the adjacencies, so that the + /// appropriate amount of memory can be set aside for storage when + /// the adjacency structure is created. + /// + /// The triangulation is assumed to involve 3-node triangles. + /// + /// Two nodes are "adjacent" if they are both nodes in some triangle. + /// Also, a node is considered to be adjacent to itself. + /// + /// Diagram: + /// + /// 3 + /// s |\ + /// i | \ + /// d | \ + /// e | \ side 1 + /// | \ + /// 2 | \ + /// | \ + /// 1-------2 + /// + /// side 3 + /// + int[] AdjacencyCount(Mesh mesh) + { + int i; + int node; + int n1, n2, n3; + int tri_id; + int neigh_id; + + int[] adj_rows = new int[node_num + 1]; + + // Set every node to be adjacent to itself. + for (node = 0; node < node_num; node++) + { + adj_rows[node] = 1; + } + + // Examine each triangle. + foreach (var tri in mesh.triangles.Values) + { + tri_id = tri.id; + + n1 = tri.vertices[0].id; + n2 = tri.vertices[1].id; + n3 = tri.vertices[2].id; + + // Add edge (1,2) if this is the first occurrence, that is, if + // the edge (1,2) is on a boundary (nid <= 0) or if this triangle + // is the first of the pair in which the edge occurs (tid < nid). + neigh_id = tri.neighbors[2].triangle.id; + + if (neigh_id < 0 || tri_id < neigh_id) + { + adj_rows[n1] += 1; + adj_rows[n2] += 1; + } + + // Add edge (2,3). + neigh_id = tri.neighbors[0].triangle.id; + + if (neigh_id < 0 || tri_id < neigh_id) + { + adj_rows[n2] += 1; + adj_rows[n3] += 1; + } + + // Add edge (3,1). + neigh_id = tri.neighbors[1].triangle.id; + + if (neigh_id < 0 || tri_id < neigh_id) + { + adj_rows[n3] += 1; + adj_rows[n1] += 1; + } + } + + // We used ADJ_COL to count the number of entries in each column. + // Convert it to pointers into the ADJ array. + for (node = node_num; 1 <= node; node--) + { + adj_rows[node] = adj_rows[node - 1]; + } + + adj_rows[0] = 1; + for (i = 1; i <= node_num; i++) + { + adj_rows[i] = adj_rows[i - 1] + adj_rows[i]; + } + + return adj_rows; + } + + /// + /// Sets adjacencies in a triangulation. + /// + /// + /// This routine can be used to create the compressed column storage + /// for a linear triangle finite element discretization of Poisson's + /// equation in two dimensions. + /// + int[] AdjacencySet(Mesh mesh, int[] rows) + { + // Output list, stores the actual adjacency information. + int[] list; + + // Copy of the adjacency rows input. + int[] rowsCopy = new int[node_num]; + Array.Copy(rows, rowsCopy, node_num); + + int i, n = rows[node_num] - 1; + + list = new int[n]; + for (i = 0; i < n; i++) + { + list[i] = -1; + } + + // Set every node to be adjacent to itself. + for (i = 0; i < node_num; i++) + { + list[rowsCopy[i] - 1] = i; + rowsCopy[i] += 1; + } + + int n1, n2, n3; // Vertex numbers. + int tid, nid; // Triangle and neighbor id. + + // Examine each triangle. + foreach (var tri in mesh.triangles.Values) + { + tid = tri.id; + + n1 = tri.vertices[0].id; + n2 = tri.vertices[1].id; + n3 = tri.vertices[2].id; + + // Add edge (1,2) if this is the first occurrence, that is, if + // the edge (1,2) is on a boundary (nid <= 0) or if this triangle + // is the first of the pair in which the edge occurs (tid < nid). + nid = tri.neighbors[2].triangle.id; + + if (nid < 0 || tid < nid) + { + list[rowsCopy[n1] - 1] = n2; + rowsCopy[n1] += 1; + list[rowsCopy[n2] - 1] = n1; + rowsCopy[n2] += 1; + } + + // Add edge (2,3). + nid = tri.neighbors[0].triangle.id; + + if (nid < 0 || tid < nid) + { + list[rowsCopy[n2] - 1] = n3; + rowsCopy[n2] += 1; + list[rowsCopy[n3] - 1] = n2; + rowsCopy[n3] += 1; + } + + // Add edge (3,1). + nid = tri.neighbors[1].triangle.id; + + if (nid < 0 || tid < nid) + { + list[rowsCopy[n1] - 1] = n3; + rowsCopy[n1] += 1; + list[rowsCopy[n3] - 1] = n1; + rowsCopy[n3] += 1; + } + } + + int k1, k2; + + // Ascending sort the entries for each node. + for (i = 0; i < node_num; i++) + { + k1 = rows[i]; + k2 = rows[i + 1] - 1; + HeapSort(list, k1 - 1, k2 + 1 - k1); + } + + return list; + } + + #endregion + + #region Heap sort + + /// + /// Reorders an array of integers into a descending heap. + /// + /// the size of the input array. + /// an unsorted array. + /// + /// A heap is an array A with the property that, for every index J, + /// A[J] >= A[2*J+1] and A[J] >= A[2*J+2], (as long as the indices + /// 2*J+1 and 2*J+2 are legal). + /// + /// Diagram: + /// + /// A(0) + /// / \ + /// A(1) A(2) + /// / \ / \ + /// A(3) A(4) A(5) A(6) + /// / \ / \ + /// A(7) A(8) A(9) A(10) + /// + private void CreateHeap(int[] a, int offset, int size) + { + int i; + int ifree; + int key; + int m; + + // Only nodes (N/2)-1 down to 0 can be "parent" nodes. + for (i = (size / 2) - 1; 0 <= i; i--) + { + // Copy the value out of the parent node. + // Position IFREE is now "open". + key = a[offset + i]; + ifree = i; + + for (; ; ) + { + // Positions 2*IFREE + 1 and 2*IFREE + 2 are the descendants of position + // IFREE. (One or both may not exist because they equal or exceed N.) + m = 2 * ifree + 1; + + // Does the first position exist? + if (size <= m) + { + break; + } + else + { + // Does the second position exist? + if (m + 1 < size) + { + // If both positions exist, take the larger of the two values, + // and update M if necessary. + if (a[offset + m] < a[offset + m + 1]) + { + m = m + 1; + } + } + + // If the large descendant is larger than KEY, move it up, + // and update IFREE, the location of the free position, and + // consider the descendants of THIS position. + if (key < a[offset + m]) + { + a[offset + ifree] = a[offset + m]; + ifree = m; + } + else + { + break; + } + } + } + + // When you have stopped shifting items up, return the item you + // pulled out back to the heap. + a[offset + ifree] = key; + } + + return; + } + + + /// + /// ascending sorts an array of integers using heap sort. + /// + /// Number of entries in the array. + /// Array to be sorted; + private void HeapSort(int[] a, int offset, int size) + { + int n1; + int temp; + + if (size <= 1) + { + return; + } + + // 1: Put A into descending heap form. + CreateHeap(a, offset, size); + + // 2: Sort A. + // The largest object in the heap is in A[0]. + // Move it to position A[N-1]. + temp = a[offset]; + a[offset] = a[offset + size - 1]; + a[offset + size - 1] = temp; + + // Consider the diminished heap of size N1. + for (n1 = size - 1; 2 <= n1; n1--) + { + // Restore the heap structure of the initial N1 entries of A. + CreateHeap(a, offset, n1); + + // Take the largest object from A[0] and move it to A[N1-1]. + temp = a[offset]; + a[offset] = a[offset + n1 - 1]; + a[offset + n1 - 1] = temp; + } + + return; + } + + #endregion + } +} diff --git a/Triangle.NET/Triangle/Tools/CuthillMcKee.cs b/Triangle.NET/Triangle/Tools/CuthillMcKee.cs index 069a4a4..d5ef2da 100644 --- a/Triangle.NET/Triangle/Tools/CuthillMcKee.cs +++ b/Triangle.NET/Triangle/Tools/CuthillMcKee.cs @@ -26,16 +26,8 @@ namespace TriangleNet.Tools // Number of nodes in the mesh. int node_num; - // Pointers into the actual adjacency structure adj. Information about row k is - // stored in entries adj_row(k) through adj_row(k+1)-1 of adj. Size: node_num + 1 - int[] adj_row; - - // Number of adjacency entries. - int adj_num; - - // The adjacency structure. For each row, it contains the column indices - // of the nonzero entries. Size: adj_num - int[] adj; + // The adjacency matrix of the mesh. + AdjacencyMatrix matrix; /// /// Gets the permutation vector for the Reverse Cuthill-McKee numbering. @@ -52,13 +44,9 @@ namespace TriangleNet.Tools mesh.Renumber(NodeNumbering.Linear); // Set up the adj_row adjacency pointer array. - this.adj_row = AdjacencyCount(mesh); - this.adj_num = adj_row[node_num] - 1; + matrix = new AdjacencyMatrix(mesh); - // Set up the adj adjacency array. - this.adj = AdjacencySet(mesh, this.adj_row); - - bandwidth1 = Bandwidth(); + bandwidth1 = matrix.Bandwidth(); // Compute the RCM permutation. int[] perm = GenerateRcm(); @@ -69,237 +57,13 @@ namespace TriangleNet.Tools if (Behavior.Verbose) { - SimpleLog.Instance.Info(String.Format("Reverse Cuthill-McKee (Banwidth: {0} > {1})", + SimpleLog.Instance.Info(String.Format("Reverse Cuthill-McKee (Bandwidth: {0} > {1})", bandwidth1, bandwidth2)); } return perm_inv; } - #region Adjacency structure - - /// - /// Counts adjacencies in a triangulation. - /// - /// - /// This routine is called to count the adjacencies, so that the - /// appropriate amount of memory can be set aside for storage when - /// the adjacency structure is created. - /// - /// The triangulation is assumed to involve 3-node triangles. - /// - /// Two nodes are "adjacent" if they are both nodes in some triangle. - /// Also, a node is considered to be adjacent to itself. - /// - /// Diagram: - /// - /// 3 - /// s |\ - /// i | \ - /// d | \ - /// e | \ side 1 - /// | \ - /// 2 | \ - /// | \ - /// 1-------2 - /// - /// side 3 - /// - int[] AdjacencyCount(Mesh mesh) - { - int i; - int node; - int n1, n2, n3; - int tri_id; - int neigh_id; - - int[] adj_rows = new int[node_num + 1]; - - // Set every node to be adjacent to itself. - for (node = 0; node < node_num; node++) - { - adj_rows[node] = 1; - } - - // Examine each triangle. - foreach (var tri in mesh.triangles.Values) - { - tri_id = tri.id; - - n1 = tri.vertices[0].id; - n2 = tri.vertices[1].id; - n3 = tri.vertices[2].id; - - // Add edge (1,2) if this is the first occurrence, that is, if - // the edge (1,2) is on a boundary (nid <= 0) or if this triangle - // is the first of the pair in which the edge occurs (tid < nid). - neigh_id = tri.neighbors[2].triangle.id; - - if (neigh_id < 0 || tri_id < neigh_id) - { - adj_rows[n1] += 1; - adj_rows[n2] += 1; - } - - // Add edge (2,3). - neigh_id = tri.neighbors[0].triangle.id; - - if (neigh_id < 0 || tri_id < neigh_id) - { - adj_rows[n2] += 1; - adj_rows[n3] += 1; - } - - // Add edge (3,1). - neigh_id = tri.neighbors[1].triangle.id; - - if (neigh_id < 0 || tri_id < neigh_id) - { - adj_rows[n3] += 1; - adj_rows[n1] += 1; - } - } - - // We used ADJ_COL to count the number of entries in each column. - // Convert it to pointers into the ADJ array. - for (node = node_num; 1 <= node; node--) - { - adj_rows[node] = adj_rows[node - 1]; - } - - adj_rows[0] = 1; - for (i = 1; i <= node_num; i++) - { - adj_rows[i] = adj_rows[i - 1] + adj_rows[i]; - } - - return adj_rows; - } - - /// - /// Sets adjacencies in a triangulation. - /// - /// - /// This routine can be used to create the compressed column storage - /// for a linear triangle finite element discretization of Poisson's - /// equation in two dimensions. - /// - int[] AdjacencySet(Mesh mesh, int[] rows) - { - // Output list, stores the actual adjacency information. - int[] list; - - // Copy of the adjacency rows input. - int[] rowsCopy = new int[node_num]; - Array.Copy(rows, rowsCopy, node_num); - - int i, n = rows[node_num] - 1; - - list = new int[n]; - for (i = 0; i < n; i++) - { - list[i] = -1; - } - - // Set every node to be adjacent to itself. - for (i = 0; i < node_num; i++) - { - list[rowsCopy[i] - 1] = i; - rowsCopy[i] += 1; - } - - int n1, n2, n3; // Vertex numbers. - int tid, nid; // Triangle and neighbor id. - - // Examine each triangle. - foreach (var tri in mesh.triangles.Values) - { - tid = tri.id; - - n1 = tri.vertices[0].id; - n2 = tri.vertices[1].id; - n3 = tri.vertices[2].id; - - // Add edge (1,2) if this is the first occurrence, that is, if - // the edge (1,2) is on a boundary (nid <= 0) or if this triangle - // is the first of the pair in which the edge occurs (tid < nid). - nid = tri.neighbors[2].triangle.id; - - if (nid < 0 || tid < nid) - { - list[rowsCopy[n1] - 1] = n2; - rowsCopy[n1] += 1; - list[rowsCopy[n2] - 1] = n1; - rowsCopy[n2] += 1; - } - - // Add edge (2,3). - nid = tri.neighbors[0].triangle.id; - - if (nid < 0 || tid < nid) - { - list[rowsCopy[n2] - 1] = n3; - rowsCopy[n2] += 1; - list[rowsCopy[n3] - 1] = n2; - rowsCopy[n3] += 1; - } - - // Add edge (3,1). - nid = tri.neighbors[1].triangle.id; - - if (nid < 0 || tid < nid) - { - list[rowsCopy[n1] - 1] = n3; - rowsCopy[n1] += 1; - list[rowsCopy[n3] - 1] = n1; - rowsCopy[n3] += 1; - } - } - - int k1, k2; - - // Ascending sort the entries for each node. - for (i = 0; i < node_num; i++) - { - k1 = rows[i]; - k2 = rows[i + 1] - 1; - HeapSort(list, k1 - 1, k2 + 1 - k1); - } - - return list; - } - - /// - /// Computes the bandwidth of an adjacency matrix. - /// - /// Bandwidth of the adjacency matrix. - int Bandwidth() - { - int band_hi; - int band_lo; - int col; - int i; - int j; - int value; - - band_lo = 0; - band_hi = 0; - - for (i = 0; i < node_num; i++) - { - for (j = adj_row[i]; j <= adj_row[i + 1] - 1; j++) - { - col = adj[j - 1]; - band_lo = Math.Max(band_lo, i - col); - band_hi = Math.Max(band_hi, col - i); - } - } - - value = band_lo + 1 + band_hi; - - return value; - } - /// /// Computes the bandwidth of a permuted adjacency matrix. /// @@ -312,7 +76,9 @@ namespace TriangleNet.Tools /// int PermBandwidth(int[] perm, int[] perm_inv) { - int bandwidth; + int[] adj_row = matrix.AdjacencyRow; + int[] adj = matrix.Adjacency; + int col, i, j; int band_lo = 0; @@ -328,12 +94,8 @@ namespace TriangleNet.Tools } } - bandwidth = band_lo + 1 + band_hi; - - return bandwidth; + return band_lo + 1 + band_hi; } - - #endregion #region RCM @@ -421,6 +183,9 @@ namespace TriangleNet.Tools /// void Rcm(int root, int[] mask, int[] perm, int offset, ref int iccsze) { + int[] adj_row = matrix.AdjacencyRow; + int[] adj = matrix.Adjacency; + int fnbr; int i, j, k, l; int jstop, jstrt; @@ -572,6 +337,9 @@ namespace TriangleNet.Tools void FindRoot(ref int root, int[] mask, ref int level_num, int[] level_row, int[] level, int offset) { + int[] adj_row = matrix.AdjacencyRow; + int[] adj = matrix.Adjacency; + int iccsze; int j, jstrt; int k, kstop, kstrt; @@ -682,6 +450,9 @@ namespace TriangleNet.Tools void GetLevelSet(ref int root, int[] mask, ref int level_num, int[] level_row, int[] level, int offset) { + int[] adj_row = matrix.AdjacencyRow; + int[] adj = matrix.Adjacency; + int i, iccsze; int j, jstop, jstrt; int lbegin, lvlend, lvsize; @@ -766,6 +537,9 @@ namespace TriangleNet.Tools /// void Degree(int root, int[] mask, int[] deg, ref int iccsze, int[] ls, int offset) { + int[] adj_row = matrix.AdjacencyRow; + int[] adj = matrix.Adjacency; + int i, ideg; int j, jstop, jstrt; int lbegin, lvlend; @@ -880,129 +654,6 @@ namespace TriangleNet.Tools return; } - /// - /// Reorders an array of integers into a descending heap. - /// - /// the size of the input array. - /// an unsorted array. - /// - /// A heap is an array A with the property that, for every index J, - /// A[J] >= A[2*J+1] and A[J] >= A[2*J+2], (as long as the indices - /// 2*J+1 and 2*J+2 are legal). - /// - /// Diagram: - /// - /// A(0) - /// / \ - /// A(1) A(2) - /// / \ / \ - /// A(3) A(4) A(5) A(6) - /// / \ / \ - /// A(7) A(8) A(9) A(10) - /// - void CreateHeap(int[] a, int offset, int size) - { - int i; - int ifree; - int key; - int m; - - // Only nodes (N/2)-1 down to 0 can be "parent" nodes. - for (i = (size / 2) - 1; 0 <= i; i--) - { - // Copy the value out of the parent node. - // Position IFREE is now "open". - key = a[offset + i]; - ifree = i; - - for (; ; ) - { - // Positions 2*IFREE + 1 and 2*IFREE + 2 are the descendants of position - // IFREE. (One or both may not exist because they equal or exceed N.) - m = 2 * ifree + 1; - - // Does the first position exist? - if (size <= m) - { - break; - } - else - { - // Does the second position exist? - if (m + 1 < size) - { - // If both positions exist, take the larger of the two values, - // and update M if necessary. - if (a[offset + m] < a[offset + m + 1]) - { - m = m + 1; - } - } - - // If the large descendant is larger than KEY, move it up, - // and update IFREE, the location of the free position, and - // consider the descendants of THIS position. - if (key < a[offset + m]) - { - a[offset + ifree] = a[offset + m]; - ifree = m; - } - else - { - break; - } - } - } - - // When you have stopped shifting items up, return the item you - // pulled out back to the heap. - a[offset + ifree] = key; - } - - return; - } - - - /// - /// ascending sorts an array of integers using heap sort. - /// - /// Number of entries in the array. - /// Array to be sorted; - void HeapSort(int[] a, int offset, int size) - { - int n1; - int temp; - - if (size <= 1) - { - return; - } - - // 1: Put A into descending heap form. - CreateHeap(a, offset, size); - - // 2: Sort A. - // The largest object in the heap is in A[0]. - // Move it to position A[N-1]. - temp = a[offset]; - a[offset] = a[offset + size - 1]; - a[offset + size - 1] = temp; - - // Consider the diminished heap of size N1. - for (n1 = size - 1; 2 <= n1; n1--) - { - // Restore the heap structure of the initial N1 entries of A. - CreateHeap(a, offset, n1); - - // Take the largest object from A[0] and move it to A[N1-1]. - temp = a[offset]; - a[offset] = a[offset + n1 - 1]; - a[offset + n1 - 1] = temp; - } - - return; - } - #endregion } } diff --git a/Triangle.NET/Triangle/Triangle.csproj b/Triangle.NET/Triangle/Triangle.csproj index 35c3610..0cdc54f 100644 --- a/Triangle.NET/Triangle/Triangle.csproj +++ b/Triangle.NET/Triangle/Triangle.csproj @@ -82,6 +82,7 @@ +