diff --git a/Triangle.NET/TestApp/Controls/AngleHistogram.cs b/Triangle.NET/TestApp/Controls/AngleHistogram.cs index 8bc8135..3837e40 100644 --- a/Triangle.NET/TestApp/Controls/AngleHistogram.cs +++ b/Triangle.NET/TestApp/Controls/AngleHistogram.cs @@ -95,6 +95,7 @@ namespace MeshExplorer.Controls if (maxAngleCount == 0) { this.maxAngles = null; + return; } this.Invalidate(); diff --git a/Triangle.NET/TestApp/Controls/DarkCheckBox.cs b/Triangle.NET/TestApp/Controls/DarkCheckBox.cs index eb786c9..1932145 100644 --- a/Triangle.NET/TestApp/Controls/DarkCheckBox.cs +++ b/Triangle.NET/TestApp/Controls/DarkCheckBox.cs @@ -82,6 +82,37 @@ namespace MeshExplorer.Controls InitializeComponent(); } + private void DrawText(Graphics g, Color forecolor, Point location) + { + if (this.UseCompatibleTextRendering) + { + //using (StringFormat stringFormat = this.CreateStringFormat()) + { + if (this.Enabled) + { + g.DrawString(this.Text, base.Font, new SolidBrush(forecolor), location.X, location.Y); + } + else + { + g.DrawString(this.Text, base.Font, new SolidBrush(forecolor), location.X, location.Y); + } + } + } + else + { + //TextFormatFlags textFormatFlags = this.CreateTextFormatFlags(); + if (this.Enabled) + { + TextRenderer.DrawText((IDeviceContext)g, this.Text, this.Font, location, forecolor); + } + else + { + //forecolor = TextRenderer.DisabledTextColor(this.BackColor); + TextRenderer.DrawText((IDeviceContext)g, this.Text, this.Font, location, forecolor); + } + } + } + #region Control overrides protected override void OnPaint(PaintEventArgs e) @@ -153,7 +184,7 @@ namespace MeshExplorer.Controls e.Graphics.TextRenderingHint = TextRenderingHint.ClearTypeGridFit; SizeF szL = e.Graphics.MeasureString(this.Text, base.Font, this.Width); - e.Graphics.DrawString(this.Text, base.Font, new SolidBrush(text_color), boxSize + 4, (this.Height - szL.Height) / 2); + DrawText(e.Graphics, text_color, new Point(boxSize + 4, (int)((this.Height - szL.Height) / 2) + 1)); if (this.isChecked) { diff --git a/Triangle.NET/TestApp/Controls/DarkTabControl.cs b/Triangle.NET/TestApp/Controls/DarkTabControl.cs index addea6c..831ee07 100644 --- a/Triangle.NET/TestApp/Controls/DarkTabControl.cs +++ b/Triangle.NET/TestApp/Controls/DarkTabControl.cs @@ -16,7 +16,7 @@ namespace MeshExplorer.Controls /// /// Summary description for FlatTabControl. /// - public class DarkTabControl : System.Windows.Forms.TabControl + public class DarkTabControl : System.Windows.Forms.TabControl { #region Designer diff --git a/Triangle.NET/TestApp/FormMain.Designer.cs b/Triangle.NET/TestApp/FormMain.Designer.cs index ba37d72..74ae0f7 100644 --- a/Triangle.NET/TestApp/FormMain.Designer.cs +++ b/Triangle.NET/TestApp/FormMain.Designer.cs @@ -101,6 +101,8 @@ this.menuTools = new System.Windows.Forms.ToolStripMenuItem(); this.menuToolsGen = new System.Windows.Forms.ToolStripMenuItem(); this.menuToolsCheck = new System.Windows.Forms.ToolStripMenuItem(); + this.toolStripSeparator4 = new System.Windows.Forms.ToolStripSeparator(); + this.menuToolsRcm = new System.Windows.Forms.ToolStripMenuItem(); this.btnMesh = new MeshExplorer.Controls.DarkButton(); this.renderControl1 = new MeshExplorer.Controls.RendererControl(); ((System.ComponentModel.ISupportInitialize)(this.splitContainer1)).BeginInit(); @@ -755,7 +757,7 @@ this.label19.Name = "label19"; this.label19.Size = new System.Drawing.Size(134, 40); this.label19.TabIndex = 0; - this.label19.Text = "Beta 2 (2012-06-20)\r\nChristian Woltering\r\nMIT"; + this.label19.Text = "Beta 2 (2012-07-30)\r\nChristian Woltering\r\nMIT"; // // label18 // @@ -894,7 +896,9 @@ // this.menuTools.DropDownItems.AddRange(new System.Windows.Forms.ToolStripItem[] { this.menuToolsGen, - this.menuToolsCheck}); + this.menuToolsCheck, + this.toolStripSeparator4, + this.menuToolsRcm}); this.menuTools.Name = "menuTools"; this.menuTools.Size = new System.Drawing.Size(46, 24); this.menuTools.Text = "Tools"; @@ -902,7 +906,7 @@ // menuToolsGen // this.menuToolsGen.Name = "menuToolsGen"; - this.menuToolsGen.Size = new System.Drawing.Size(157, 22); + this.menuToolsGen.Size = new System.Drawing.Size(195, 22); this.menuToolsGen.Text = "Input Generator"; this.menuToolsGen.Click += new System.EventHandler(this.menuToolsGenerator_Click); // @@ -910,10 +914,22 @@ // this.menuToolsCheck.Enabled = false; this.menuToolsCheck.Name = "menuToolsCheck"; - this.menuToolsCheck.Size = new System.Drawing.Size(157, 22); + this.menuToolsCheck.Size = new System.Drawing.Size(195, 22); this.menuToolsCheck.Text = "Check Mesh"; this.menuToolsCheck.Click += new System.EventHandler(this.menuToolsCheck_Click); // + // toolStripSeparator4 + // + this.toolStripSeparator4.Name = "toolStripSeparator4"; + this.toolStripSeparator4.Size = new System.Drawing.Size(192, 6); + // + // menuToolsRcm + // + this.menuToolsRcm.Name = "menuToolsRcm"; + this.menuToolsRcm.Size = new System.Drawing.Size(195, 22); + this.menuToolsRcm.Text = "Renumber nodes (RCM)"; + this.menuToolsRcm.Click += new System.EventHandler(this.menuToolsRcm_Click); + // // btnMesh // this.btnMesh.Enabled = false; @@ -1050,6 +1066,8 @@ private System.Windows.Forms.Label lbNumVert2; private System.Windows.Forms.Label label23; private Controls.DarkCheckBox cbConformDel; + private System.Windows.Forms.ToolStripSeparator toolStripSeparator4; + private System.Windows.Forms.ToolStripMenuItem menuToolsRcm; } } diff --git a/Triangle.NET/TestApp/FormMain.cs b/Triangle.NET/TestApp/FormMain.cs index 89ad7c1..e328b66 100644 --- a/Triangle.NET/TestApp/FormMain.cs +++ b/Triangle.NET/TestApp/FormMain.cs @@ -387,8 +387,6 @@ namespace MeshExplorer if (sfd.ShowDialog() == DialogResult.OK) { - mesh.Renumber(); - FileProcessor.Save(sfd.FileName, mesh); } } @@ -523,6 +521,19 @@ namespace MeshExplorer UpdateLog(); } + private void Renumber() + { + if (mesh == null || settings.ExceptionThrown) return; + + bool tmp = Behavior.Verbose; + Behavior.Verbose = true; + + mesh.Renumber(NodeNumbering.CuthillMcKee); + ShowLog(); + + Behavior.Verbose = tmp; + } + private void Smooth() { if (mesh == null || settings.ExceptionThrown) return; @@ -674,5 +685,10 @@ namespace MeshExplorer { this.Close(); } + + private void menuToolsRcm_Click(object sender, EventArgs e) + { + Renumber(); + } } } diff --git a/Triangle.NET/TestApp/IO/EpsImage.cs b/Triangle.NET/TestApp/IO/EpsImage.cs index 851e334..a2e319a 100644 --- a/Triangle.NET/TestApp/IO/EpsImage.cs +++ b/Triangle.NET/TestApp/IO/EpsImage.cs @@ -261,7 +261,7 @@ namespace MeshExplorer.IO private void DrawPoints(StreamWriter eps, Mesh mesh, bool label) { - int n = mesh.NumberOfVertices; + int n = mesh.Vertices.Count; int circle_size = 1; diff --git a/Triangle.NET/TestApp/IO/FileProcessor.cs b/Triangle.NET/TestApp/IO/FileProcessor.cs index d1e03b5..28e4ee5 100644 --- a/Triangle.NET/TestApp/IO/FileProcessor.cs +++ b/Triangle.NET/TestApp/IO/FileProcessor.cs @@ -97,10 +97,6 @@ namespace MeshExplorer.IO { provider = new JsonFile(); } - else if (ext == ".dat") - { - provider = new DatFile(); - } if (provider == null) { diff --git a/Triangle.NET/TestApp/IO/Formats/DatFile.cs b/Triangle.NET/TestApp/IO/Formats/DatFile.cs deleted file mode 100644 index 095084a..0000000 --- a/Triangle.NET/TestApp/IO/Formats/DatFile.cs +++ /dev/null @@ -1,80 +0,0 @@ -// ----------------------------------------------------------------------- -// -// Christian Woltering, Triangle.NET, http://triangle.codeplex.com/ -// -// ----------------------------------------------------------------------- - -namespace MeshExplorer.IO.Formats -{ - using System; - using System.Collections.Generic; - using System.Linq; - using System.Text; - using TriangleNet.IO; - using System.IO; - using MeshExplorer.Rendering; - using TriangleNet.Geometry; - using TriangleNet; - - /// - /// Read a polygon from DAT file (used by poly2tri). - /// - public class DatFile : IMeshFile - { - /// - /// Gets the supported file extensions. - /// - public string[] Extensions - { - get { return new string[] { ".dat" }; } - } - - public bool ContainsMeshData(string filename) - { - return false; - } - - public InputGeometry Read(string filename) - { - InputGeometry data = new InputGeometry(); - - string line; - string[] split; - - using (TextReader reader = new StreamReader(filename)) - { - while ((line = reader.ReadLine()) != null) - { - split = line.Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries); - - if (split.Length == 2) - { - data.AddPoint( - double.Parse(split[0], Util.Nfi), - double.Parse(split[1], Util.Nfi)); - } - } - } - - int n = data.Count; - - for (int i = 0; i < n; i++) - { - data.AddSegment(i, (i + 1) % n); - - } - - return data; - } - - public Mesh Import(string filename) - { - throw new NotImplementedException(); - } - - public void Write(Mesh mesh, string filename) - { - throw new NotImplementedException(); - } - } -} diff --git a/Triangle.NET/TestApp/IO/Formats/JsonFile.cs b/Triangle.NET/TestApp/IO/Formats/JsonFile.cs index ade67e5..3419453 100644 --- a/Triangle.NET/TestApp/IO/Formats/JsonFile.cs +++ b/Triangle.NET/TestApp/IO/Formats/JsonFile.cs @@ -16,6 +16,7 @@ namespace MeshExplorer.IO.Formats using TriangleNet.Geometry; using TriangleNet; using System.Collections; + using TriangleNet.Data; /// /// Read and write JSON files. @@ -101,10 +102,10 @@ namespace MeshExplorer.IO.Formats { using (StreamWriter writer = new StreamWriter(filename)) { - int nv = mesh.NumberOfVertices; - int ns = mesh.NumberOfSegments; + int nv = mesh.Vertices.Count; + int ns = mesh.Segments.Count; int nh = mesh.Holes.Count; - int ne = mesh.NumberOfTriangles; + int ne = mesh.Triangles.Count; writer.Write("{"); @@ -115,6 +116,11 @@ namespace MeshExplorer.IO.Formats writer.Write("\"dim\":2"); writer.Write("}"); + if (mesh.CurrentNumbering == NodeNumbering.None) + { + mesh.Renumber(NodeNumbering.Linear); + } + // Write the coordinates if (nv > 0) { @@ -397,17 +403,47 @@ namespace MeshExplorer.IO.Formats #region Write helpers - private void WritePoints(Mesh data, StreamWriter writer, int nv) + private void WritePoints(Mesh mesh, StreamWriter writer, int nv) { - int i = 0; - - StringBuilder markers = new StringBuilder(); bool useMarkers = false; - string seperator; + StringBuilder markers; writer.Write("\"points\":{\"data\":["); - foreach (var item in data.Vertices) + + if (mesh.CurrentNumbering == NodeNumbering.Linear) + { + markers = WritePoints(writer, mesh.Vertices, nv, useMarkers); + } + else + { + Vertex[] nodes = new Vertex[mesh.Vertices.Count]; + + foreach (var node in mesh.Vertices) + { + nodes[node.ID] = node; + } + + markers = WritePoints(writer, nodes, nv, useMarkers); + } + + writer.Write("]"); + if (useMarkers) + { + writer.Write(",\"markers\":[" + markers.ToString() + "]"); + } + + // TODO: writer.Write(",\"attributes\":[]"); + writer.Write("}"); + } + + private static StringBuilder WritePoints(StreamWriter writer, IEnumerable nodes, int nv, bool useMarkers) + { + StringBuilder markers = new StringBuilder(); + + int i = 0; + string seperator; + foreach (var item in nodes) { seperator = (i == nv - 1) ? String.Empty : ", "; @@ -424,13 +460,8 @@ namespace MeshExplorer.IO.Formats i++; } - writer.Write("]"); - if (useMarkers) - { - writer.Write(",\"markers\":[" + markers.ToString() + "]"); - } - //writer.Write(",\"attributes\":[]"); - writer.Write("}"); + + return markers; } private void WriteHoles(Mesh data, StreamWriter writer, int nh) diff --git a/Triangle.NET/TestApp/IO/Formats/TriangleFile.cs b/Triangle.NET/TestApp/IO/Formats/TriangleFile.cs index 1c44c98..4ffbfda 100644 --- a/Triangle.NET/TestApp/IO/Formats/TriangleFile.cs +++ b/Triangle.NET/TestApp/IO/Formats/TriangleFile.cs @@ -57,7 +57,7 @@ namespace MeshExplorer.IO.Formats public void Write(Mesh mesh, string filename) { - if (mesh.NumberOfVertices > 0) + if (mesh.Vertices.Count > 0) { format.Write(mesh, filename); } diff --git a/Triangle.NET/TestApp/IO/RasterImage.cs b/Triangle.NET/TestApp/IO/RasterImage.cs index 005e3ce..141aa56 100644 --- a/Triangle.NET/TestApp/IO/RasterImage.cs +++ b/Triangle.NET/TestApp/IO/RasterImage.cs @@ -51,13 +51,13 @@ namespace MeshExplorer.IO Bitmap bitmap; // Check if the specified width is reasonable - if (width < 2 * Math.Sqrt(mesh.NumberOfVertices)) + if (width < 2 * Math.Sqrt(mesh.Vertices.Count)) { bitmap = new Bitmap(400, 200); Graphics g = Graphics.FromImage(bitmap); g.Clear(colors.Background); - string message = String.Format("Sorry, I won't render {0} points on such a small image!", mesh.NumberOfVertices); + string message = String.Format("Sorry, I won't render {0} points on such a small image!", mesh.Vertices.Count); SizeF sz = g.MeasureString(message, SystemFonts.DefaultFont); diff --git a/Triangle.NET/TestApp/IO/SvgImage.cs b/Triangle.NET/TestApp/IO/SvgImage.cs index 1ae0778..fec887f 100644 --- a/Triangle.NET/TestApp/IO/SvgImage.cs +++ b/Triangle.NET/TestApp/IO/SvgImage.cs @@ -153,7 +153,7 @@ namespace MeshExplorer.IO private void DrawPoints(StreamWriter svg, Mesh mesh, bool label) { - int n = mesh.NumberOfVertices; + int n = mesh.Vertices.Count; int circle_size = 1; diff --git a/Triangle.NET/TestApp/Mesh Explorer.csproj b/Triangle.NET/TestApp/Mesh Explorer.csproj index 155b194..76e463e 100644 --- a/Triangle.NET/TestApp/Mesh Explorer.csproj +++ b/Triangle.NET/TestApp/Mesh Explorer.csproj @@ -107,7 +107,6 @@ - diff --git a/Triangle.NET/TestApp/Rendering/RenderData.cs b/Triangle.NET/TestApp/Rendering/RenderData.cs index 9e2fabc..d435402 100644 --- a/Triangle.NET/TestApp/Rendering/RenderData.cs +++ b/Triangle.NET/TestApp/Rendering/RenderData.cs @@ -73,7 +73,7 @@ namespace MeshExplorer.Rendering this.Segments = null; - int n = mesh.NumberOfVertices; + int n = mesh.Vertices.Count; // Convert points to float this.Points = new PointF[n]; @@ -85,7 +85,7 @@ namespace MeshExplorer.Rendering { var segs = mesh.Segments; - List segList = new List(mesh.NumberOfSegments); + List segList = new List(mesh.Segments.Count); foreach (var seg in segs) { diff --git a/Triangle.NET/TestApp/Settings.cs b/Triangle.NET/TestApp/Settings.cs index ba9bb5d..df8d601 100644 --- a/Triangle.NET/TestApp/Settings.cs +++ b/Triangle.NET/TestApp/Settings.cs @@ -59,7 +59,7 @@ namespace MeshExplorer SfdFilterIndex = 1; OfdFilter = SfdFilter; - OfdFilter += "|Polygon data (*.dat)|*.dat"; + //OfdFilter += "|Polygon data (*.dat)|*.dat"; //OfdFilter += "|COMSOL mesh (*.mphtxt)|*.mphtxt"; //OfdFilter += "|AVS UCD data (*.ucd)|*.ucd"; //OfdFilter += "|VTK data (*.vtk)|*.vtk"; diff --git a/Triangle.NET/Triangle/BadTriQueue.cs b/Triangle.NET/Triangle/BadTriQueue.cs index 6966a81..f8c797d 100644 --- a/Triangle.NET/Triangle/BadTriQueue.cs +++ b/Triangle.NET/Triangle/BadTriQueue.cs @@ -162,10 +162,6 @@ namespace TriangleNet newbad.triangorg = enqorg; newbad.triangdest = enqdest; - Vertex org = enqtri.Org(); - Vertex dest = enqtri.Dest(); - Vertex apex = enqtri.Apex(); - Enqueue(newbad); } diff --git a/Triangle.NET/Triangle/Behavior.cs b/Triangle.NET/Triangle/Behavior.cs index ea76342..38b8404 100644 --- a/Triangle.NET/Triangle/Behavior.cs +++ b/Triangle.NET/Triangle/Behavior.cs @@ -12,7 +12,7 @@ namespace TriangleNet /// /// Controls the behavior of the meshing software. /// - class Behavior + public class Behavior { /// /// Load behavior defaults. diff --git a/Triangle.NET/Triangle/Data/Vertex.cs b/Triangle.NET/Triangle/Data/Vertex.cs index 9394318..4b25ffd 100644 --- a/Triangle.NET/Triangle/Data/Vertex.cs +++ b/Triangle.NET/Triangle/Data/Vertex.cs @@ -27,18 +27,51 @@ namespace TriangleNet.Data internal VertexType type; internal Otri tri; + /// + /// Initializes a new instance of the class. + /// public Vertex() - : this(0, 0, 0) - { } + : this(0, 0, 0, 0) + { + } + /// + /// Initializes a new instance of the class. + /// + /// The x coordinate of the vertex. + /// The y coordinate of the vertex. public Vertex(double x, double y) - : this(x, y, 0) - { } + : this(x, y, 0, 0) + { + } + /// + /// Initializes a new instance of the class. + /// + /// The x coordinate of the vertex. + /// The y coordinate of the vertex. + /// The boundary mark. public Vertex(double x, double y, int mark) + : this(x, y, mark, 0) + { + } + + /// + /// Initializes a new instance of the class. + /// + /// The x coordinate of the vertex. + /// The y coordinate of the vertex. + /// The boundary mark. + /// The number of point attributes. + public Vertex(double x, double y, int mark, int attribs) : base(x, y, mark) { this.type = VertexType.InputVertex; + + if (attribs > 0) + { + this.attributes = new double[attribs]; + } } #region Public properties diff --git a/Triangle.NET/Triangle/Enums.cs b/Triangle.NET/Triangle/Enums.cs index 786b7b1..e0d227c 100644 --- a/Triangle.NET/Triangle/Enums.cs +++ b/Triangle.NET/Triangle/Enums.cs @@ -95,5 +95,5 @@ namespace TriangleNet /// /// Node renumbering algorithms. /// - public enum NodeNumbering { Linear, CuthillMcKee }; + public enum NodeNumbering { None, Linear, CuthillMcKee }; } diff --git a/Triangle.NET/Triangle/Geometry/InputGeometry.cs b/Triangle.NET/Triangle/Geometry/InputGeometry.cs index 7d44ec8..66afa94 100644 --- a/Triangle.NET/Triangle/Geometry/InputGeometry.cs +++ b/Triangle.NET/Triangle/Geometry/InputGeometry.cs @@ -23,6 +23,9 @@ namespace TriangleNet.Geometry BoundingBox bounds; + // Used to check consitent use of point attributes. + private int pointAttributes = -1; + /// /// Initializes a new instance of the class. /// @@ -44,6 +47,8 @@ namespace TriangleNet.Geometry regions = new List(); bounds = new BoundingBox(); + + pointAttributes = -1; } /// @@ -55,7 +60,7 @@ namespace TriangleNet.Geometry } /// - /// Gets a value indicating whether the geometry should be treated as a PLSG. + /// Indicates, whether the geometry should be treated as a PSLG. /// public bool HasSegments { @@ -95,7 +100,7 @@ namespace TriangleNet.Geometry } /// - /// Gets the list of input holes. + /// Gets the list of regions. /// public IEnumerable Regions { @@ -111,6 +116,8 @@ namespace TriangleNet.Geometry segments.Clear(); holes.Clear(); regions.Clear(); + + pointAttributes = -1; } /// @@ -136,6 +143,45 @@ namespace TriangleNet.Geometry bounds.Update(x, y); } + /// + /// Adds a point to the geometry. + /// + /// X coordinate. + /// Y coordinate. + /// Boundary marker. + /// Point attribute. + public void AddPoint(double x, double y, int boundary, double attribute) + { + AddPoint(x, y, 0, new double[] { attribute }); + } + + /// + /// Adds a point to the geometry. + /// + /// X coordinate. + /// Y coordinate. + /// Boundary marker. + /// Point attributes. + public void AddPoint(double x, double y, int boundary, double[] attribs) + { + if (pointAttributes < 0) + { + pointAttributes = attribs == null ? 0 : attribs.Length; + } + else if (attribs == null && pointAttributes > 0) + { + throw new ArgumentException("Inconsitent use of point attributes."); + } + else if (attribs != null && pointAttributes != attribs.Length) + { + throw new ArgumentException("Inconsitent use of point attributes."); + } + + points.Add(new Vertex(x, y, boundary) { attributes = attribs }); + + bounds.Update(x, y); + } + /// /// Adds a hole location to the geometry. /// diff --git a/Triangle.NET/Triangle/Geometry/Point.cs b/Triangle.NET/Triangle/Geometry/Point.cs index ce2bc0b..8dabcc2 100644 --- a/Triangle.NET/Triangle/Geometry/Point.cs +++ b/Triangle.NET/Triangle/Geometry/Point.cs @@ -1,6 +1,5 @@ // ----------------------------------------------------------------------- // -// Original Triangle code by Jonathan Richard Shewchuk, http://www.cs.cmu.edu/~quake/triangle.html // Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/ // // ----------------------------------------------------------------------- @@ -144,6 +143,11 @@ namespace TriangleNet.Geometry return (x < other.x || (x == other.x && y < other.y)) ? -1 : 1; } + public override int GetHashCode() + { + return x.GetHashCode() ^ y.GetHashCode(); + } + public override string ToString() { return String.Format("[{0},{1}]", x, y); diff --git a/Triangle.NET/Triangle/Geometry/RegionPointer.cs b/Triangle.NET/Triangle/Geometry/RegionPointer.cs index 9e730cc..3267ac8 100644 --- a/Triangle.NET/Triangle/Geometry/RegionPointer.cs +++ b/Triangle.NET/Triangle/Geometry/RegionPointer.cs @@ -1,6 +1,5 @@ // ----------------------------------------------------------------------- // -// Original Triangle code by Jonathan Richard Shewchuk, http://www.cs.cmu.edu/~quake/triangle.html // Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/ // // ----------------------------------------------------------------------- diff --git a/Triangle.NET/Triangle/IO/DataReader.cs b/Triangle.NET/Triangle/IO/DataReader.cs index 3586820..b0b0d32 100644 --- a/Triangle.NET/Triangle/IO/DataReader.cs +++ b/Triangle.NET/Triangle/IO/DataReader.cs @@ -92,7 +92,7 @@ namespace TriangleNet.IO // Create the subsegments. for (i = 0; i < mesh.insegments; i++) { - mesh.MakeSubseg(ref subseg); + mesh.MakeSegment(ref subseg); // Mark the subsegment as living. //subseg.ss.subsegs[0].ss = subseg.ss; } diff --git a/Triangle.NET/Triangle/IO/FileReader.cs b/Triangle.NET/Triangle/IO/FileReader.cs index 07deddf..254a283 100644 --- a/Triangle.NET/Triangle/IO/FileReader.cs +++ b/Triangle.NET/Triangle/IO/FileReader.cs @@ -57,30 +57,31 @@ namespace TriangleNet.IO /// The input geometry. /// The current vertex index. /// The current line. - /// Number of point attributes - static void ReadVertex(InputGeometry data, int index, string[] line, int n) + /// Number of point attributes + /// Number of point markers (0 or 1) + static void ReadVertex(InputGeometry data, int index, string[] line, int attributes, int marks) { double x = double.Parse(line[1], nfi); double y = double.Parse(line[2], nfi); int mark = 0; + double[] attribs = attributes == 0 ? null : new double[attributes]; // Read the vertex attributes. - for (int j = 0; j < n; j++) + for (int j = 0; j < attributes; j++) { if (line.Length > 3 + j) { - // TODO: - //vertex.attributes[j] = double.Parse(line[3 + j]); + attribs[j] = double.Parse(line[3 + j]); } } // Read a vertex marker. - if (line.Length > 3 + n) + if (marks > 0 && line.Length > 3 + attributes) { - mark = int.Parse(line[3 + n]); + mark = int.Parse(line[3 + attributes]); } - data.AddPoint(x, y, mark); + data.AddPoint(x, y, mark, attribs); } #endregion @@ -218,7 +219,7 @@ namespace TriangleNet.IO startIndex = int.Parse(line[0], nfi); } - ReadVertex(data, i, line, attributes); + ReadVertex(data, i, line, attributes, nodemarkers); } } } @@ -325,7 +326,7 @@ namespace TriangleNet.IO startIndex = int.Parse(line[0], nfi); } - ReadVertex(data, i, line, attributes); + ReadVertex(data, i, line, attributes, nodemarkers); } } else diff --git a/Triangle.NET/Triangle/IO/FileWriter.cs b/Triangle.NET/Triangle/IO/FileWriter.cs index a9dd932..7abcd58 100644 --- a/Triangle.NET/Triangle/IO/FileWriter.cs +++ b/Triangle.NET/Triangle/IO/FileWriter.cs @@ -12,6 +12,7 @@ namespace TriangleNet.IO using System.Globalization; using TriangleNet.Data; using TriangleNet.Geometry; + using System.Collections.Generic; /// /// Helper methods for writing Triangle file formats. @@ -25,12 +26,10 @@ namespace TriangleNet.IO /// /// /// - public static void WriteNodes(Mesh mesh, string filename) + public static void Write(Mesh mesh, string filename) { - using (StreamWriter writer = new StreamWriter(filename)) - { - FileWriter.WriteNodes(mesh, writer); - } + FileWriter.WritePoly(mesh, Path.ChangeExtension(filename, ".poly")); + FileWriter.WriteElements(mesh, Path.ChangeExtension(filename, ".ele")); } /// @@ -38,10 +37,20 @@ namespace TriangleNet.IO /// /// /// - private static void WriteNodes(Mesh mesh, StreamWriter writer) + public static void WriteNodes(Mesh mesh, string filename) { - Vertex vertex; - long outvertices = mesh.vertices.Count; + using (StreamWriter writer = new StreamWriter(filename)) + { + FileWriter.WriteNodes(writer, mesh); + } + } + + /// + /// Number the vertices and write them to a .node file. + /// + private static void WriteNodes(StreamWriter writer, Mesh mesh) + { + int outvertices = mesh.vertices.Count; Behavior behavior = mesh.behavior; @@ -50,8 +59,6 @@ namespace TriangleNet.IO outvertices = mesh.vertices.Count - mesh.undeads; } - int index = 0; - if (writer != null) { // Number of vertices, number of dimensions, number of vertex attributes, @@ -59,32 +66,69 @@ namespace TriangleNet.IO writer.WriteLine("{0} {1} {2} {3}", outvertices, mesh.mesh_dim, mesh.nextras, behavior.UseBoundaryMarkers ? "1" : "0"); - foreach (var item in mesh.vertices.Values) + if (mesh.numbering == NodeNumbering.None) { - vertex = item; + // If the mesh isn't numbered yet, use linear node numbering. + mesh.Renumber(); + } - if (!behavior.Jettison || vertex.type != VertexType.UndeadVertex) + if (mesh.numbering == NodeNumbering.Linear) + { + // If numbering is linear, just use the dictionary values. + WriteNodes(writer, mesh.vertices.Values, behavior.UseBoundaryMarkers, + mesh.nextras, behavior.Jettison); + } + else + { + // If numbering is not linear, a simple 'foreach' traversal of the dictionary + // values doesn't reflect the actual numbering. Use an array instead. + + // TODO: Could use a custom sorting function on dictionary values instead. + Vertex[] nodes = new Vertex[mesh.vertices.Count]; + + foreach (var node in mesh.vertices.Values) { - // Vertex number, x and y coordinates. - writer.Write("{0} {1} {2}", index, vertex.x.ToString(nfi), vertex.y.ToString(nfi)); - - // Write attributes. - for (int j = 0; j < mesh.nextras; j++) - { - writer.Write(" {0}", vertex.attributes[j].ToString(nfi)); - } - - if (behavior.UseBoundaryMarkers) - { - // Write the boundary marker. - writer.Write(" {0}", vertex.mark); - } - - writer.WriteLine(); - - // Assign array index to vertex ID for later use. - vertex.id = index++; + nodes[node.id] = node; } + + WriteNodes(writer, nodes, behavior.UseBoundaryMarkers, + mesh.nextras, behavior.Jettison); + } + } + } + + /// + /// Write the vertices to a stream. + /// + /// + /// + private static void WriteNodes(StreamWriter writer, IEnumerable nodes, bool markers, + int attribs, bool jettison) + { + int index = 0; + + foreach (var vertex in nodes) + { + if (!jettison || vertex.type != VertexType.UndeadVertex) + { + // Vertex number, x and y coordinates. + writer.Write("{0} {1} {2}", index, vertex.x.ToString(nfi), vertex.y.ToString(nfi)); + + // Write attributes. + for (int j = 0; j < attribs; j++) + { + writer.Write(" {0}", vertex.attributes[j].ToString(nfi)); + } + + if (markers) + { + // Write the boundary marker. + writer.Write(" {0}", vertex.mark); + } + + writer.WriteLine(); + + index++; } } } @@ -163,7 +207,7 @@ namespace TriangleNet.IO if (writeNodes) { // Write nodes to this file. - FileWriter.WriteNodes(mesh, writer); + FileWriter.WriteNodes(writer, mesh); } else { @@ -202,10 +246,11 @@ namespace TriangleNet.IO } // Holes + j = 0; writer.WriteLine("{0}", mesh.holes.Count); foreach (var hole in mesh.holes) { - writer.WriteLine("{0} {1}", hole.X.ToString(nfi), hole.Y.ToString(nfi)); + writer.WriteLine("{0} {1} {2}", j++, hole.X.ToString(nfi), hole.Y.ToString(nfi)); } // Regions diff --git a/Triangle.NET/Triangle/Mesh.cs b/Triangle.NET/Triangle/Mesh.cs index 814f0cb..9f68495 100644 --- a/Triangle.NET/Triangle/Mesh.cs +++ b/Triangle.NET/Triangle/Mesh.cs @@ -32,7 +32,6 @@ namespace TriangleNet // Stack that maintains a list of recently flipped triangles. Stack flipstack; - //FlipStacker lastflip; // TODO: Check if custom hashmap implementation could be faster. @@ -69,12 +68,11 @@ namespace TriangleNet // Triangular bounding box vertices. internal Vertex infvertex1, infvertex2, infvertex3; - // The 'triangle' that occupies all of "outer space." + // The 'triangle' that occupies all of 'outer space'. internal static Triangle dummytri; - // The omnipresent subsegment. Referenced by any triangle or - // subsegment that isn't really connected to a subsegment at - // that location. + // The omnipresent subsegment. Referenced by any triangle or subsegment + // that isn't really connected to a subsegment at that location. internal static Segment dummysub; // Pointer to a recently visited triangle. Improves point location if @@ -84,6 +82,9 @@ namespace TriangleNet // Controls the behavior of the mesh instance. internal Behavior behavior; + // The current node numbering + internal NodeNumbering numbering; + #endregion #region Public properties @@ -99,7 +100,7 @@ namespace TriangleNet /// /// Gets the mesh vertices. /// - public IEnumerable Vertices + public ICollection Vertices { get { return this.vertices.Values; } } @@ -115,7 +116,7 @@ namespace TriangleNet /// /// Gets the mesh triangles. /// - public IEnumerable Triangles + public ICollection Triangles { get { return this.triangles.Values; } } @@ -123,7 +124,7 @@ namespace TriangleNet /// /// Gets the mesh segments. /// - public IEnumerable Segments + public ICollection Segments { get { return this.subsegs.Values; } } @@ -133,21 +134,6 @@ namespace TriangleNet /// public int NumberOfInputPoints { get { return invertices; } } - /// - /// Gets the number of mesh vertices. - /// - public int NumberOfVertices { get { return this.vertices.Count; } } - - /// - /// Gets the number of mesh triangles. - /// - public int NumberOfTriangles { get { return this.triangles.Count; } } - - /// - /// Gets the number of mesh segments. - /// - public int NumberOfSegments { get { return this.subsegs.Count; } } - /// /// Gets the number of mesh edges. /// @@ -158,6 +144,14 @@ namespace TriangleNet /// public bool IsPolygon { get { return this.insegments > 0; } } + /// + /// Gets the current node numbering. + /// + public NodeNumbering CurrentNumbering + { + get { return numbering; } + } + #endregion /// @@ -461,6 +455,8 @@ namespace TriangleNet /// public void Smooth() { + numbering = NodeNumbering.None; + //ISmoother smoother = new CvdSmoother(this); //smoother.Smooth(); } @@ -478,29 +474,38 @@ namespace TriangleNet /// public void Renumber(NodeNumbering num) { + // Don't need to do anything if the nodes are already numbered. + if (num == this.numbering) + { + return; + } + int id; if (num == NodeNumbering.Linear) { id = 0; - foreach (var item in this.vertices.Values) + foreach (var node in this.vertices.Values) { - item.id = id++; + node.id = id++; } } - else + else if (num == NodeNumbering.CuthillMcKee) { - //CuthillMcKee rcm = new CuthillMcKee(); - //int[] perm_inv = rcm.Renumber(this); + CuthillMcKee rcm = new CuthillMcKee(); + int[] perm_inv = rcm.Renumber(this); - //// Permute the node indices. - //foreach (var node in this.vertices.Values) - //{ - // node.id = perm_inv[node.id]; - //} + // Permute the node indices. + foreach (var node in this.vertices.Values) + { + node.id = perm_inv[node.id]; + } } - // Triangles will always be numbered from 0..n + // Remember the current numbering. + numbering = num; + + // Triangles will always be numbered from 0 to n-1 id = 0; foreach (var item in this.triangles.Values) { @@ -723,6 +728,8 @@ namespace TriangleNet /// private void Reset() { + numbering = NodeNumbering.None; + recenttri.triangle = null; // No triangle has been visited yet. undeads = 0; // No eliminated input vertices yet. checksegments = false; // There are no segments in the triangulation yet. @@ -816,13 +823,10 @@ namespace TriangleNet throw new Exception("Input must have at least three input vertices."); } - this.nextras = 0; // TODO: points[0].Attributes == null ? 0 : points[0].Attributes.Length; + this.nextras = points[0].attributes == null ? 0 : points[0].attributes.Length; foreach (Vertex vertex in points) { - // TODO: Set vertex attributes. - //vertex.attribs = points[i].Attributes; - vertex.hash = this.hash_vtx++; vertex.id = vertex.hash; @@ -856,7 +860,7 @@ namespace TriangleNet /// Create a new subsegment with orientation zero. /// /// Reference to the new subseg. - internal void MakeSubseg(ref Osub newsubseg) + internal void MakeSegment(ref Osub newsubseg) { Segment seg = new Segment(); seg.hash = this.hash_seg++; @@ -866,6 +870,7 @@ namespace TriangleNet subsegs.Add(seg.hash, seg); } + #endregion #region Manipulation @@ -1485,7 +1490,7 @@ namespace TriangleNet if (newsubseg.seg == dummysub) { // Make new subsegment and initialize its vertices. - MakeSubseg(ref newsubseg); + MakeSegment(ref newsubseg); newsubseg.SetOrg(tridest); newsubseg.SetDest(triorg); newsubseg.SetSegOrg(tridest); @@ -2510,7 +2515,7 @@ namespace TriangleNet Vertex leftvertex, rightvertex; Vertex newvertex; InsertVertexResult success; - FindDirectionResult collinear; + double ex, ey; double tx, ty; double etx, ety; @@ -2535,13 +2540,23 @@ namespace TriangleNet throw new Exception("Attempt to find intersection of parallel segments."); } split = (ey * etx - ex * ety) / denom; + // Create the new vertex. - newvertex = new Vertex(torg.x + split * (tdest.x - torg.x), - torg.y + split * (tdest.y - torg.y), splitsubseg.seg.boundary); + newvertex = new Vertex( + torg.x + split * (tdest.x - torg.x), + torg.y + split * (tdest.y - torg.y), + splitsubseg.seg.boundary, + this.nextras); + newvertex.hash = this.hash_vtx++; newvertex.id = newvertex.hash; - // TODO: nextras //(vertex) poolalloc(&m.vertices); + // Interpolate its attributes. + for (int i = 0; i < nextras; i++) + { + newvertex.attributes[i] = torg.attributes[i] + split * (tdest.attributes[i] - torg.attributes[i]); + } + vertices.Add(newvertex.hash, newvertex); // Insert the intersection vertex. This should always succeed. @@ -2576,7 +2591,8 @@ namespace TriangleNet // Inserting the vertex may have caused edge flips. We wish to rediscover // the edge connecting endpoint1 to the new intersection vertex. - collinear = FindDirection(ref splittri, endpoint1); + FindDirection(ref splittri, endpoint1); + rightvertex = splittri.Dest(); leftvertex = splittri.Apex(); if ((leftvertex.x == endpoint1.x) && (leftvertex.y == endpoint1.y)) diff --git a/Triangle.NET/Triangle/NewLocation.cs b/Triangle.NET/Triangle/NewLocation.cs index 62213d4..a07fc53 100644 --- a/Triangle.NET/Triangle/NewLocation.cs +++ b/Triangle.NET/Triangle/NewLocation.cs @@ -34,28 +34,25 @@ namespace TriangleNet /// /// Find a new location for a Steiner point. /// - /// /// /// /// - /// /// /// /// /// - public void FindLocation(Vertex torg, Vertex tdest, Vertex tapex, - Vertex circumcenter, ref double xi, ref double eta, bool offcenter, Otri badotri) // TODO: ref circumcenter??? + /// + public Point FindLocation(Vertex torg, Vertex tdest, Vertex tapex, + ref double xi, ref double eta, bool offcenter, Otri badotri) { // Based on using -U switch, call the corresponding function if (behavior.MaxAngle == 0.0) { - FindNewLocationWithoutMaxAngle(torg, tdest, tapex, circumcenter, ref xi, ref eta, true, badotri); - } - else - { - FindNewLocation(torg, tdest, tapex, circumcenter, ref xi, ref eta, true, badotri); + return FindNewLocationWithoutMaxAngle(torg, tdest, tapex, ref xi, ref eta, true, badotri); } + // With max angle + return FindNewLocation(torg, tdest, tapex, ref xi, ref eta, true, badotri); } /// @@ -69,8 +66,8 @@ namespace TriangleNet /// /// /// - private void FindNewLocationWithoutMaxAngle(Vertex torg, Vertex tdest, Vertex tapex, - Vertex circumcenter, ref double xi, ref double eta, bool offcenter, Otri badotri) + private Point FindNewLocationWithoutMaxAngle(Vertex torg, Vertex tdest, Vertex tapex, + ref double xi, ref double eta, bool offcenter, Otri badotri) { double offconstant = behavior.Offconstant; @@ -740,6 +737,8 @@ namespace TriangleNet }// end of relocation }// end of almostGood + Point circumcenter = new Point(); + if (relocated <= 0) { circumcenter.x = torg.x + dx; @@ -754,6 +753,7 @@ namespace TriangleNet xi = (yao * dx - xao * dy) * (2.0 * denominator); eta = (xdo * dy - ydo * dx) * (2.0 * denominator); + return circumcenter; } /// @@ -767,8 +767,8 @@ namespace TriangleNet /// /// /// - private void FindNewLocation(Vertex torg, Vertex tdest, Vertex tapex, - Vertex circumcenter, ref double xi, ref double eta, bool offcenter, Otri badotri) + private Point FindNewLocation(Vertex torg, Vertex tdest, Vertex tapex, + ref double xi, ref double eta, bool offcenter, Otri badotri) { double offconstant = behavior.Offconstant; @@ -1300,9 +1300,7 @@ namespace TriangleNet && (IsBadTriangleAngle(middleAngleCorner.x, middleAngleCorner.y, largestAngleCorner.x, largestAngleCorner.y, petal_slab_inter_x_first, petal_slab_inter_y_first)) && MinDistanceToNeighbor(petal_slab_inter_x_first, petal_slab_inter_y_first, ref neighborotri) > MinDistanceToNeighbor(line_inter_x, line_inter_y, ref neighborotri)) { - // - /// check the neighbor's vertices also, which one if better - + // check the neighbor's vertices also, which one if better //slab and petal intersection is advised dxFirstSuggestion = petal_slab_inter_x_first - torg.x; dyFirstSuggestion = petal_slab_inter_y_first - torg.y; @@ -1329,26 +1327,19 @@ namespace TriangleNet // go back to circumcenter dxFirstSuggestion = dx; dyFirstSuggestion = dy; - - } else { // intersection point is suggested dxFirstSuggestion = line_inter_x - torg.x; dyFirstSuggestion = line_inter_y - torg.y; - - } - - } else {// we are not creating a bad triangle // slab intersection is advised dxFirstSuggestion = line_result[2] - torg.x; dyFirstSuggestion = line_result[3] - torg.y; - } } //------------------------------------------------------// @@ -1361,8 +1352,6 @@ namespace TriangleNet // go back to circumcenter dxFirstSuggestion = dx; dyFirstSuggestion = dy; - - } else { @@ -1370,11 +1359,8 @@ namespace TriangleNet // neighbor's circumcenter is suggested dxFirstSuggestion = voronoiOrInter[2] - torg.x; dyFirstSuggestion = voronoiOrInter[3] - torg.y; - - } } - } else { // there is no voronoi vertex between intersection point and circumcenter @@ -1397,7 +1383,6 @@ namespace TriangleNet //slab and petal intersection is advised dxFirstSuggestion = petal_slab_inter_x_first - torg.x; dyFirstSuggestion = petal_slab_inter_y_first - torg.y; - } else { // slab intersection point is further away @@ -1421,32 +1406,25 @@ namespace TriangleNet // go back to circumcenter dxFirstSuggestion = dx; dyFirstSuggestion = dy; - } else { // intersection point is suggested dxFirstSuggestion = line_inter_x - torg.x; dyFirstSuggestion = line_inter_y - torg.y; - } - - } else {// we are not creating a bad triangle // slab intersection is advised dxFirstSuggestion = line_result[2] - torg.x; dyFirstSuggestion = line_result[3] - torg.y; - } } //------------------------------------------------------// - } else { - if (IsBadTriangleAngle(largestAngleCorner.x, largestAngleCorner.y, middleAngleCorner.x, middleAngleCorner.y, inter_x, inter_y)) { //printf("testtriangle returned false! bad triangle\n"); @@ -1469,16 +1447,12 @@ namespace TriangleNet // go back to circumcenter dxFirstSuggestion = dx; dyFirstSuggestion = dy; - - - } else { // intersection point is suggested dxFirstSuggestion = inter_x - torg.x; dyFirstSuggestion = inter_y - torg.y; - } } else @@ -1486,7 +1460,6 @@ namespace TriangleNet // intersection point is suggested dxFirstSuggestion = inter_x - torg.x; dyFirstSuggestion = inter_y - torg.y; - } } } @@ -1614,8 +1587,6 @@ namespace TriangleNet // slab and petal intersection is advised dxSecondSuggestion = petal_slab_inter_x_second - torg.x; dySecondSuggestion = petal_slab_inter_y_second - torg.y; - - } else { // slab intersection point is further away @@ -1639,8 +1610,6 @@ namespace TriangleNet // go back to circumcenter dxSecondSuggestion = dx; dySecondSuggestion = dy; - - } else { @@ -1649,15 +1618,12 @@ namespace TriangleNet dySecondSuggestion = line_inter_y - torg.y; } - - } else {// we are not creating a bad triangle // slab intersection is advised dxSecondSuggestion = line_result[2] - torg.x; dySecondSuggestion = line_result[3] - torg.y; - } } //------------------------------------------------------// @@ -1669,16 +1635,12 @@ namespace TriangleNet // go back to circumcenter dxSecondSuggestion = dx; dySecondSuggestion = dy; - - } else { // we are not creating a bad triangle // neighbor's circumcenter is suggested dxSecondSuggestion = voronoiOrInter[2] - torg.x; dySecondSuggestion = voronoiOrInter[3] - torg.y; - - } } } @@ -1703,7 +1665,6 @@ namespace TriangleNet // slab and petal intersection is advised dxSecondSuggestion = petal_slab_inter_x_second - torg.x; dySecondSuggestion = petal_slab_inter_y_second - torg.y; - } else { // slab intersection point is further away ; @@ -1727,28 +1688,23 @@ namespace TriangleNet // go back to circumcenter dxSecondSuggestion = dx; dySecondSuggestion = dy; - } else { // intersection point is suggested dxSecondSuggestion = line_inter_x - torg.x; dySecondSuggestion = line_inter_y - torg.y; - } - - } else - {// we are not creating a bad triangle + { + // we are not creating a bad triangle // slab intersection is advised dxSecondSuggestion = line_result[2] - torg.x; dySecondSuggestion = line_result[3] - torg.y; - } } //------------------------------------------------------// - } else { @@ -1773,24 +1729,19 @@ namespace TriangleNet // go back to circumcenter dxSecondSuggestion = dx; dySecondSuggestion = dy; - - } else { // intersection point is suggested dxSecondSuggestion = inter_x - torg.x; dySecondSuggestion = inter_y - torg.y; - } } else { - // intersection point is suggested dxSecondSuggestion = inter_x - torg.x; dySecondSuggestion = inter_y - torg.y; - } } } @@ -1827,13 +1778,11 @@ namespace TriangleNet { dx = dxSecondSuggestion; dy = dySecondSuggestion; - } else { dx = dxFirstSuggestion; dy = dyFirstSuggestion; - } } else if (neighborNotFound_first) @@ -1850,13 +1799,11 @@ namespace TriangleNet { dx = dxSecondSuggestion; dy = dySecondSuggestion; - } else { dx = dxFirstSuggestion; dy = dyFirstSuggestion; - } } else if (neighborNotFound_second) @@ -1873,13 +1820,11 @@ namespace TriangleNet { dx = dxSecondSuggestion; dy = dySecondSuggestion; - } else { dx = dxFirstSuggestion; dy = dyFirstSuggestion; - } } else @@ -1896,16 +1841,13 @@ namespace TriangleNet { dx = dxSecondSuggestion; dy = dySecondSuggestion; - } else { dx = dxFirstSuggestion; dy = dyFirstSuggestion; - } } - } else { // acute : consider other direction @@ -1923,13 +1865,11 @@ namespace TriangleNet { dx = dxSecondSuggestion; dy = dySecondSuggestion; - } else { dx = dxFirstSuggestion; dy = dyFirstSuggestion; - } } else if (neighborNotFound_first) @@ -1946,13 +1886,11 @@ namespace TriangleNet { dx = dxSecondSuggestion; dy = dySecondSuggestion; - } else { dx = dxFirstSuggestion; dy = dyFirstSuggestion; - } } else if (neighborNotFound_second) @@ -1969,13 +1907,11 @@ namespace TriangleNet { dx = dxSecondSuggestion; dy = dySecondSuggestion; - } else { dx = dxFirstSuggestion; dy = dyFirstSuggestion; - } } else @@ -1992,13 +1928,11 @@ namespace TriangleNet { dx = dxSecondSuggestion; dy = dySecondSuggestion; - } else { dx = dxFirstSuggestion; dy = dyFirstSuggestion; - } } @@ -2006,6 +1940,8 @@ namespace TriangleNet }// end of relocation }// end of almostGood + Point circumcenter = new Point(); + if (relocated <= 0) { circumcenter.x = torg.x + dx; @@ -2019,6 +1955,7 @@ namespace TriangleNet xi = (yao * dx - xao * dy) * (2.0 * denominator); eta = (xdo * dy - ydo * dx) * (2.0 * denominator); + return circumcenter; } /// diff --git a/Triangle.NET/Triangle/Primitives.cs b/Triangle.NET/Triangle/Primitives.cs index e5b1db5..1f40079 100644 --- a/Triangle.NET/Triangle/Primitives.cs +++ b/Triangle.NET/Triangle/Primitives.cs @@ -19,10 +19,9 @@ namespace TriangleNet { static double splitter; // Used to split double factors for exact multiplication. static double epsilon; // Floating-point machine epsilon. - static double resulterrbound; - static double ccwerrboundA, ccwerrboundB, ccwerrboundC; - static double iccerrboundA, iccerrboundB, iccerrboundC; - static double o3derrboundA, o3derrboundB, o3derrboundC; + //static double resulterrbound; + static double ccwerrboundA; // ccwerrboundB, ccwerrboundC; + static double iccerrboundA; // iccerrboundB, iccerrboundC; /// /// Initialize the variables used for exact arithmetic. @@ -69,16 +68,13 @@ namespace TriangleNet } while ((check != 1.0) && (check != lastcheck)); splitter += 1.0; // Error bounds for orientation and incircle tests. - resulterrbound = (3.0 + 8.0 * epsilon) * epsilon; + //resulterrbound = (3.0 + 8.0 * epsilon) * epsilon; ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon; - ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon; - ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon; + //ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon; + //ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon; iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon; - iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon; - iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon; - o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon; - o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon; - o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon; + //iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon; + //iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon; } /// @@ -386,7 +382,7 @@ namespace TriangleNet ref double xi, ref double eta) { double xdo, ydo, xao, yao; - double dodist, aodist, dadist; + double dodist, aodist; double denominator; double dx, dy; @@ -399,8 +395,6 @@ namespace TriangleNet yao = tapex.y - torg.y; dodist = xdo * xdo + ydo * ydo; aodist = xao * xao + yao * yao; - dadist = (tdest.x - tapex.x) * (tdest.x - tapex.x) + - (tdest.y - tapex.y) * (tdest.y - tapex.y); if (Behavior.NoExact) { diff --git a/Triangle.NET/Triangle/Quality.cs b/Triangle.NET/Triangle/Quality.cs index 1aa06f3..4ad3afe 100644 --- a/Triangle.NET/Triangle/Quality.cs +++ b/Triangle.NET/Triangle/Quality.cs @@ -373,7 +373,7 @@ namespace TriangleNet double area; double dist1, dist2; - double maxedge, maxangle; + double maxangle; torg = testtri.Org(); tdest = testtri.Dest(); @@ -463,21 +463,21 @@ namespace TriangleNet if ((apexlen > orglen) && (apexlen > destlen)) { // The edge opposite the apex is longest. - maxedge = apexlen; + // maxedge = apexlen; // Find the cosine of the angle at the apex. maxangle = (orglen + destlen - apexlen) / (2 * Math.Sqrt(orglen) * Math.Sqrt(destlen)); } else if (orglen > destlen) { // The edge opposite the origin is longest. - maxedge = orglen; + // maxedge = orglen; // Find the cosine of the angle at the origin. maxangle = (apexlen + destlen - orglen) / (2 * Math.Sqrt(apexlen) * Math.Sqrt(destlen)); } else { // The edge opposite the destination is longest. - maxedge = destlen; + // maxedge = destlen; // Find the cosine of the angle at the destination. maxangle = (apexlen + orglen - destlen) / (2 * Math.Sqrt(apexlen) * Math.Sqrt(orglen)); } @@ -600,7 +600,6 @@ namespace TriangleNet double split; double multiplier, divisor; bool acuteorg, acuteorg2, acutedest, acutedest2; - int dummy; // Note that steinerleft == -1 if an unlimited number // of Steiner points is allowed. @@ -727,7 +726,13 @@ namespace TriangleNet } // Create the new vertex. - newvertex = new Vertex(); // TODO: mesh.nextras + newvertex = new Vertex( + eorg.x + split * (edest.x - eorg.x), + eorg.y + split * (edest.y - eorg.y), + currentenc.Mark(), + mesh.nextras); + + newvertex.type = VertexType.SegmentVertex; newvertex.hash = mesh.hash_vtx++; newvertex.id = newvertex.hash; @@ -741,9 +746,6 @@ namespace TriangleNet + split * (edest.attributes[i] - eorg.attributes[i]); } - newvertex.x = eorg.x + split * (edest.x - eorg.x); - newvertex.y = eorg.y + split * (edest.y - eorg.y); - if (!Behavior.NoExact) { // Roundoff in the above calculation may yield a 'newvertex' @@ -764,9 +766,6 @@ namespace TriangleNet } } - newvertex.mark = currentenc.Mark(); - newvertex.type = VertexType.SegmentVertex; - // Check whether the new vertex lies on an endpoint. if (((newvertex.x == eorg.x) && (newvertex.y == eorg.y)) || ((newvertex.x == edest.x) && (newvertex.y == edest.y))) @@ -791,9 +790,9 @@ namespace TriangleNet mesh.steinerleft--; } // Check the two new subsegments to see if they're encroached. - dummy = CheckSeg4Encroach(ref currentenc); + CheckSeg4Encroach(ref currentenc); currentenc.NextSelf(); - dummy = CheckSeg4Encroach(ref currentenc); + CheckSeg4Encroach(ref currentenc); } // Set subsegment's origin to NULL. This makes it possible to detect dead @@ -828,7 +827,7 @@ namespace TriangleNet { Otri badotri = default(Otri); Vertex borg, bdest, bapex; - Vertex newvertex; + Point newloc; // Location of the new vertex double xi = 0, eta = 0; InsertVertexResult success; bool errorflag; @@ -846,7 +845,6 @@ namespace TriangleNet { errorflag = false; // Create a new vertex at the triangle's circumcenter. - newvertex = new Vertex(); // TODO: mesh.nextras // Using the original (simpler) Steiner point location method // for mesh refinement. @@ -854,19 +852,17 @@ namespace TriangleNet // reset VertexType? if (behavior.FixedArea || behavior.VarArea) { - Point tmp = Primitives.FindCircumcenter(borg, bdest, bapex, ref xi, ref eta, behavior.Offconstant); - newvertex.x = tmp.x; - newvertex.y = tmp.y; + newloc = Primitives.FindCircumcenter(borg, bdest, bapex, ref xi, ref eta, behavior.Offconstant); } else { - newLocation.FindLocation(borg, bdest, bapex, newvertex, ref xi, ref eta, true, badotri); + newloc = newLocation.FindLocation(borg, bdest, bapex, ref xi, ref eta, true, badotri); } // Check whether the new vertex lies on a triangle vertex. - if (((newvertex.x == borg.x) && (newvertex.y == borg.y)) || - ((newvertex.x == bdest.x) && (newvertex.y == bdest.y)) || - ((newvertex.x == bapex.x) && (newvertex.y == bapex.y))) + if (((newloc.x == borg.x) && (newloc.y == borg.y)) || + ((newloc.x == bdest.x) && (newloc.y == bdest.y)) || + ((newloc.x == bapex.x) && (newloc.y == bapex.y))) { if (Behavior.Verbose) { @@ -876,6 +872,11 @@ namespace TriangleNet } else { + // The new vertex must be in the interior, and therefore is a + // free vertex with a marker of zero. + Vertex newvertex = new Vertex(newloc.x, newloc.y, 0, mesh.nextras); + newvertex.type = VertexType.FreeVertex; + for (int i = 0; i < mesh.nextras; i++) { // Interpolate the vertex attributes at the circumcenter. @@ -883,10 +884,6 @@ namespace TriangleNet + xi * (bdest.attributes[i] - borg.attributes[i]) + eta * (bapex.attributes[i] - borg.attributes[i]); } - // The new vertex must be in the interior, and therefore is a - // free vertex with a marker of zero. - newvertex.type = VertexType.FreeVertex; - //newvertex.mark = 0; // Ensure that the handle 'badotri' does not represent the longest // edge of the triangle. This ensures that the circumcenter must diff --git a/Triangle.NET/Triangle/Tools/CuthillMcKee.cs b/Triangle.NET/Triangle/Tools/CuthillMcKee.cs new file mode 100644 index 0000000..069a4a4 --- /dev/null +++ b/Triangle.NET/Triangle/Tools/CuthillMcKee.cs @@ -0,0 +1,1008 @@ +// ----------------------------------------------------------------------- +// +// Original Matlab code by John Burkardt, Florida State University +// 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; + using TriangleNet.Log; + + /// + /// Applies the Cuthill and McKee renumbering algorithm to reduce the bandwidth of + /// the adjacency matrix associated with the mesh. + /// + /// + /// Some useful slides: + /// http://bobbyness.net/NerdyStuff/node%20ordering/node_ordering.html + /// + public class CuthillMcKee + { + // 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; + + /// + /// Gets the permutation vector for the Reverse Cuthill-McKee numbering. + /// + /// The mesh. + /// Permutation vector. + public int[] Renumber(Mesh mesh) + { + int bandwidth1, bandwidth2; + + this.node_num = mesh.vertices.Count; + + // Algorithm needs linear numbering of the nodes. + 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; + + // Set up the adj adjacency array. + this.adj = AdjacencySet(mesh, this.adj_row); + + bandwidth1 = Bandwidth(); + + // Compute the RCM permutation. + int[] perm = GenerateRcm(); + + int[] perm_inv = PermInverse(node_num, perm); + + bandwidth2 = PermBandwidth(perm, perm_inv); + + if (Behavior.Verbose) + { + SimpleLog.Instance.Info(String.Format("Reverse Cuthill-McKee (Banwidth: {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. + /// + /// The permutation. + /// The inverse permutation. + /// Bandwidth of the permuted adjacency matrix. + /// + /// The matrix is defined by the adjacency information and a permutation. + /// The routine also computes the bandwidth and the size of the envelope. + /// + int PermBandwidth(int[] perm, int[] perm_inv) + { + int bandwidth; + int col, i, j; + + int band_lo = 0; + int band_hi = 0; + + for (i = 0; i < node_num; i++) + { + for (j = adj_row[perm[i]]; j <= adj_row[perm[i] + 1] - 1; j++) + { + col = perm_inv[adj[j - 1]]; + band_lo = Math.Max(band_lo, i - col); + band_hi = Math.Max(band_hi, col - i); + } + } + + bandwidth = band_lo + 1 + band_hi; + + return bandwidth; + } + + #endregion + + #region RCM + + /// + /// Finds the reverse Cuthill-Mckee ordering for a general graph. + /// + /// The RCM ordering. + /// + /// For each connected component in the graph, the routine obtains + /// an ordering by calling RCM. + /// + int[] GenerateRcm() + { + int[] perm = new int[node_num]; + + int i, num, root; + int iccsze = 0; + int level_num = 0; + + /// Index vector for a level structure. The level structure is stored in the + /// currently unused spaces in the permutation vector PERM. + int[] level_row = new int[node_num + 1]; + + /// Marks variables that have been numbered. + int[] mask = new int[node_num]; + + for (i = 0; i < node_num; i++) + { + mask[i] = 1; + } + + num = 1; + + for (i = 0; i < node_num; i++) + { + // For each masked connected component... + if (mask[i] != 0) + { + root = i; + + // Find a pseudo-peripheral node ROOT. The level structure found by + // ROOT_FIND is stored starting at PERM(NUM). + FindRoot(ref root, mask, ref level_num, level_row, perm, num - 1); + + // RCM orders the component using ROOT as the starting node. + Rcm(root, mask, perm, num - 1, ref iccsze); + + num += iccsze; + + // We can stop once every node is in one of the connected components. + if (node_num < num) + { + return perm; + } + } + } + + return perm; + } + + /// + /// RCM renumbers a connected component by the reverse Cuthill McKee algorithm. + /// + /// the node that defines the connected component. It is used as the starting + /// point for the RCM ordering. + /// Input/output, int MASK(NODE_NUM), a mask for the nodes. Only those nodes with + /// nonzero input mask values are considered by the routine. The nodes numbered by RCM will have + /// their mask values set to zero. + /// Output, int PERM(NODE_NUM), the RCM ordering. + /// Output, int ICCSZE, the size of the connected component that has been numbered. + /// the number of nodes. + /// + /// The connected component is specified by a node ROOT and a mask. + /// The numbering starts at the root node. + /// + /// An outline of the algorithm is as follows: + /// + /// X(1) = ROOT. + /// + /// for ( I = 1 to N-1) + /// Find all unlabeled neighbors of X(I), + /// assign them the next available labels, in order of increasing degree. + /// + /// When done, reverse the ordering. + /// + void Rcm(int root, int[] mask, int[] perm, int offset, ref int iccsze) + { + int fnbr; + int i, j, k, l; + int jstop, jstrt; + int lbegin, lnbr, lperm, lvlend; + int nbr, node; + + /// 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. + int[] deg = new int[node_num]; + + // Find the degrees of the nodes in the component specified by MASK and ROOT. + Degree(root, mask, deg, ref iccsze, perm, offset); + + mask[root] = 0; + + if (iccsze <= 1) + { + return; + } + + lvlend = 0; + lnbr = 1; + + // LBEGIN and LVLEND point to the beginning and + // the end of the current level respectively. + while (lvlend < lnbr) + { + lbegin = lvlend + 1; + lvlend = lnbr; + + for (i = lbegin; i <= lvlend; i++) + { + // For each node in the current level... + node = perm[offset + i - 1]; + jstrt = adj_row[node]; + jstop = adj_row[node + 1] - 1; + + // Find the unnumbered neighbors of NODE. + + // FNBR and LNBR point to the first and last neighbors + // of the current node in PERM. + fnbr = lnbr + 1; + + for (j = jstrt; j <= jstop; j++) + { + nbr = adj[j - 1]; + + if (mask[nbr] != 0) + { + lnbr += 1; + mask[nbr] = 0; + perm[offset + lnbr - 1] = nbr; + } + } + + // Node has neighbors + if (lnbr > fnbr) + { + // Sort the neighbors of NODE in increasing order by degree. + // Linear insertion is used. + k = fnbr; + + while (k < lnbr) + { + l = k; + k = k + 1; + nbr = perm[offset + k - 1]; + + while (fnbr < l) + { + lperm = perm[offset + l - 1]; + + if (deg[lperm - 1] <= deg[nbr - 1]) + { + break; + } + + perm[offset + l] = lperm; + l = l - 1; + } + perm[offset + l] = nbr; + } + } + } + } + + // We now have the Cuthill-McKee ordering. Reverse it. + ReverseVector(perm, offset, iccsze); + + return; + } + + /// + /// Finds a pseudo-peripheral node. + /// + /// On input, ROOT is a node in the the component of the graph for + /// which a pseudo-peripheral node is sought. On output, ROOT is the pseudo-peripheral + /// node obtained. + /// MASK[NODE_NUM], specifies a section subgraph. Nodes for which MASK + /// is zero are ignored by FNROOT. + /// Output, int LEVEL_NUM, is the number of levels in the level + /// structure rooted at the node ROOT. + /// Output, int LEVEL_ROW(NODE_NUM+1), the level structure array pair + /// containing the level structure found. + /// Output, int LEVEL(NODE_NUM), the level structure array pair + /// containing the level structure found. + /// the number of nodes. + /// + /// The diameter of a graph is the maximum distance (number of edges) + /// between any two nodes of the graph. + /// + /// The eccentricity of a node is the maximum distance between that + /// node and any other node of the graph. + /// + /// A peripheral node is a node whose eccentricity equals the + /// diameter of the graph. + /// + /// A pseudo-peripheral node is an approximation to a peripheral node; + /// it may be a peripheral node, but all we know is that we tried our + /// best. + /// + /// The routine is given a graph, and seeks pseudo-peripheral nodes, + /// using a modified version of the scheme of Gibbs, Poole and + /// Stockmeyer. It determines such a node for the section subgraph + /// specified by MASK and ROOT. + /// + /// The routine also determines the level structure associated with + /// the given pseudo-peripheral node; that is, how far each node + /// is from the pseudo-peripheral node. The level structure is + /// returned as a list of nodes LS, and pointers to the beginning + /// of the list of nodes that are at a distance of 0, 1, 2, ..., + /// NODE_NUM-1 from the pseudo-peripheral node. + /// + /// Reference: + /// Alan George, Joseph Liu, + /// Computer Solution of Large Sparse Positive Definite Systems, + /// Prentice Hall, 1981. + /// + /// Norman Gibbs, William Poole, Paul Stockmeyer, + /// An Algorithm for Reducing the Bandwidth and Profile of a Sparse Matrix, + /// SIAM Journal on Numerical Analysis, + /// Volume 13, pages 236-250, 1976. + /// + /// Norman Gibbs, + /// Algorithm 509: A Hybrid Profile Reduction Algorithm, + /// ACM Transactions on Mathematical Software, + /// Volume 2, pages 378-387, 1976. + /// + void FindRoot(ref int root, int[] mask, ref int level_num, int[] level_row, + int[] level, int offset) + { + int iccsze; + int j, jstrt; + int k, kstop, kstrt; + int mindeg; + int nghbor, ndeg; + int node; + int level_num2 = 0; + + // Determine the level structure rooted at ROOT. + GetLevelSet(ref root, mask, ref level_num, level_row, level, offset); + + // Count the number of nodes in this level structure. + iccsze = level_row[level_num] - 1; + + // Extreme cases: + // A complete graph has a level set of only a single level. + // Every node is equally good (or bad). + // or + // A "line graph" 0--0--0--0--0 has every node in its only level. + // By chance, we've stumbled on the ideal root. + if (level_num == 1 || level_num == iccsze) + { + return; + } + + // Pick any node from the last level that has minimum degree + // as the starting point to generate a new level set. + for (; ; ) + { + mindeg = iccsze; + + jstrt = level_row[level_num - 1]; + root = level[offset + jstrt - 1]; + + if (jstrt < iccsze) + { + for (j = jstrt; j <= iccsze; j++) + { + node = level[offset + j - 1]; + ndeg = 0; + kstrt = adj_row[node - 1]; + kstop = adj_row[node] - 1; + + for (k = kstrt; k <= kstop; k++) + { + nghbor = adj[k - 1]; + if (mask[nghbor] > 0) + { + ndeg += 1; + } + } + + if (ndeg < mindeg) + { + root = node; + mindeg = ndeg; + } + } + } + + // Generate the rooted level structure associated with this node. + GetLevelSet(ref root, mask, ref level_num2, level_row, level, offset); + + // If the number of levels did not increase, accept the new ROOT. + if (level_num2 <= level_num) + { + break; + } + + level_num = level_num2; + + // In the unlikely case that ROOT is one endpoint of a line graph, + // we can exit now. + if (iccsze <= level_num) + { + break; + } + } + + return; + } + + /// + /// Generates the connected level structure rooted at a given node. + /// + /// the node at which the level structure is to be rooted. + /// MASK[NODE_NUM]. On input, only nodes with nonzero MASK are to be processed. + /// On output, those nodes which were included in the level set have MASK set to 1. + /// Output, int LEVEL_NUM, the number of levels in the level structure. ROOT is + /// in level 1. The neighbors of ROOT are in level 2, and so on. + /// Output, int LEVEL_ROW[NODE_NUM+1], the rooted level structure. + /// Output, int LEVEL[NODE_NUM], the rooted level structure. + /// the number of nodes. + /// + /// Only nodes for which MASK is nonzero will be considered. + /// + /// The root node chosen by the user is assigned level 1, and masked. + /// All (unmasked) nodes reachable from a node in level 1 are + /// assigned level 2 and masked. The process continues until there + /// are no unmasked nodes adjacent to any node in the current level. + /// The number of levels may vary between 2 and NODE_NUM. + /// + /// Reference: + /// Alan George, Joseph Liu, + /// Computer Solution of Large Sparse Positive Definite Systems, + /// Prentice Hall, 1981. + /// + void GetLevelSet(ref int root, int[] mask, ref int level_num, int[] level_row, + int[] level, int offset) + { + int i, iccsze; + int j, jstop, jstrt; + int lbegin, lvlend, lvsize; + int nbr; + int node; + + mask[root] = 0; + level[offset] = root; + level_num = 0; + lvlend = 0; + iccsze = 1; + + // LBEGIN is the pointer to the beginning of the current level, and + // LVLEND points to the end of this level. + for (; ; ) + { + lbegin = lvlend + 1; + lvlend = iccsze; + level_num += 1; + level_row[level_num - 1] = lbegin; + + // Generate the next level by finding all the masked neighbors of nodes + // in the current level. + for (i = lbegin; i <= lvlend; i++) + { + node = level[offset + i - 1]; + jstrt = adj_row[node]; + jstop = adj_row[node + 1] - 1; + + for (j = jstrt; j <= jstop; j++) + { + nbr = adj[j - 1]; + + if (mask[nbr] != 0) + { + iccsze += 1; + level[offset + iccsze - 1] = nbr; + mask[nbr] = 0; + } + } + } + + // Compute the current level width (the number of nodes encountered.) + // If it is positive, generate the next level. + lvsize = iccsze - lvlend; + + if (lvsize <= 0) + { + break; + } + } + + level_row[level_num] = lvlend + 1; + + // Reset MASK to 1 for the nodes in the level structure. + for (i = 0; i < iccsze; i++) + { + mask[level[offset + i]] = 1; + } + + return; + } + + /// + /// Computes the degrees of the nodes in the connected component. + /// + /// the node that defines the connected component. + /// MASK[NODE_NUM], is nonzero for those nodes which are to be considered. + /// Output, int DEG[NODE_NUM], contains, for each node in the connected component, its degree. + /// Output, int ICCSIZE, the number of nodes in the connected component. + /// Output, int LS[NODE_NUM], stores in entries 1 through ICCSIZE the nodes in the + /// connected component, starting with ROOT, and proceeding by levels. + /// the number of nodes. + /// + /// The connected component is specified by MASK and ROOT. + /// Nodes for which MASK is zero are ignored. + /// + /// Reference: + /// Alan George, Joseph Liu, + /// Computer Solution of Large Sparse Positive Definite Systems, + /// Prentice Hall, 1981. + /// + void Degree(int root, int[] mask, int[] deg, ref int iccsze, int[] ls, int offset) + { + int i, ideg; + int j, jstop, jstrt; + int lbegin, lvlend; + int lvsize = 1; + int nghbr, node; + + // The sign of ADJ_ROW(I) is used to indicate if node I has been considered. + ls[offset] = root; + adj_row[root] = -adj_row[root]; + lvlend = 0; + iccsze = 1; + + // If the current level width is nonzero, generate another level. + while (lvsize > 0) + { + // LBEGIN is the pointer to the beginning of the current level, and + // LVLEND points to the end of this level. + lbegin = lvlend + 1; + lvlend = iccsze; + + // Find the degrees of nodes in the current level, + // and at the same time, generate the next level. + for (i = lbegin; i <= lvlend; i++) + { + node = ls[offset + i - 1]; + jstrt = -adj_row[node]; + jstop = Math.Abs(adj_row[node + 1]) - 1; + ideg = 0; + + for (j = jstrt; j <= jstop; j++) + { + nghbr = adj[j - 1]; + + if (mask[nghbr] != 0) // EDIT: [nbr - 1] + { + ideg = ideg + 1; + + if (0 <= adj_row[nghbr]) // EDIT: [nbr - 1] + { + adj_row[nghbr] = -adj_row[nghbr]; // EDIT: [nbr - 1] + iccsze = iccsze + 1; + ls[offset + iccsze - 1] = nghbr; + } + } + } + deg[node] = ideg; + } + + // Compute the current level width. + lvsize = iccsze - lvlend; + } + + // Reset ADJ_ROW to its correct sign and return. + for (i = 0; i < iccsze; i++) + { + node = ls[offset + i]; + adj_row[node] = -adj_row[node]; + } + + return; + } + + #endregion + + #region Tools + + /// + /// Produces the inverse of a given permutation. + /// + /// Number of items permuted. + /// PERM[N], a permutation. + /// The inverse permutation. + int[] PermInverse(int n, int[] perm) + { + int[] perm_inv = new int[node_num]; + + int i; + + for (i = 0; i < n; i++) + { + perm_inv[perm[i]] = i; + } + + return perm_inv; + } + + /// + /// Reverses the elements of an integer vector. + /// + /// number of entries in the array. + /// the array to be reversed. + /// + /// Input: + /// N = 5, + /// A = ( 11, 12, 13, 14, 15 ). + /// + /// Output: + /// A = ( 15, 14, 13, 12, 11 ). + /// + void ReverseVector(int[] a, int offset, int size) + { + int i; + int j; + + for (i = 0; i < size / 2; i++) + { + j = a[offset + i]; + a[offset + i] = a[offset + size - 1 - i]; + a[offset + size - 1 - i] = j; + } + + 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/Tools/QualityMeasure.cs b/Triangle.NET/Triangle/Tools/QualityMeasure.cs index cffd14e..9c8acfc 100644 --- a/Triangle.NET/Triangle/Tools/QualityMeasure.cs +++ b/Triangle.NET/Triangle/Tools/QualityMeasure.cs @@ -1,5 +1,6 @@ // ----------------------------------------------------------------------- // +// Original Matlab code by John Burkardt, Florida State University // Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/ // // ----------------------------------------------------------------------- @@ -239,6 +240,8 @@ namespace TriangleNet.Tools /// public int Bandwidth() { + if (mesh == null) return 0; + // Lower and upper bandwidth of the matrix int ml = 0, mu = 0; diff --git a/Triangle.NET/Triangle/Triangle.csproj b/Triangle.NET/Triangle/Triangle.csproj index cdb523a..b0ccecb 100644 --- a/Triangle.NET/Triangle/Triangle.csproj +++ b/Triangle.NET/Triangle/Triangle.csproj @@ -82,6 +82,7 @@ +