TriangleNet.Examples Quadratic Elements

QuadraticElements.cs

namespace TriangleNet.Examples
{
    using System;
    using System.Collections.Generic;
    using System.Linq;
    using TriangleNet;
    using TriangleNet.Geometry;

    /// <summary>
    /// Create extra nodes for quadratic sub-parametric elements.
    /// </summary>
    public class QuadraticElements
    {
        /// <summary>
        /// Gets the list of new vertices (edge midpoints).
        /// </summary>
        public List<Point> Vertices { get; private set; }

        /// <summary>
        /// Gets the new vertex indices for each triangle.
        /// </summary>
        /// <remarks>
        /// For a triangle with index k, Indices[k, i] (i = 0, 1, 2) corresponds to
        /// the vertex on the edge shared with the i-th neighbor, e.g. Indices[k, 0]
        /// is the vertex shared by triangle k and it's neighbor N0 and lies on the
        /// side opposite of P0.
        /// </remarks>
        public int[,] Indices { get; private set; }

        /// <summary>
        /// Initializes a new instance of the <see cref="QuadraticElements" /> class.
        /// </summary>
        /// <param name="mesh">The mesh.</param>
        public QuadraticElements(Mesh mesh)
        {
            QuadraticOrder(mesh);
        }

        private void QuadraticOrder(Mesh mesh)
        {
            mesh.Renumber();

            Vertices = new List<Point>(mesh.NumberOfEdges);
            Indices = new int[mesh.Triangles.Count(), 3];

            Point p, q, middle;
            var pts = new Point[3]; // Triangle vertices
            ITriangle neighbor, neighbor2;
            int nid, nid2;
            int orient = 0, count = 0;

            foreach (var tri in mesh.Triangles)
            {
                pts[0] = tri.GetVertex(0);
                pts[1] = tri.GetVertex(1);
                pts[2] = tri.GetVertex(2);

                for (int i = 0; i < 3; i++)
                {
                    // The neighbor opposite of vertex i.
                    GetNeighbor(tri, i, out neighbor, out nid);

                    // Consider each edge only once.
                    if ((tri.ID < nid) || (nid < 0))
                    {
                        // Get the edge endpoints.
                        p = pts[(i + 1) % 3];
                        q = pts[(i + 2) % 3];

                        // Create a new node in the middle of the edge.
                        middle = new Point(0.5 * (p.X + q.X), 0.5 * (p.Y + q.Y),
                            Math.Min(p.Label, q.Label));

                        Vertices.Add(middle);

                        // Record the new node. 
                        Indices[tri.ID, i] = count;

                        // Record the new node in the neighbor element. 
                        if (nid != -1)
                        {
                            GetNeighbor(neighbor, orient, out neighbor2, out nid2);

                            // Get the neighbors orientation.
                            while (nid2 != tri.ID)
                            {
                                orient = (orient + 1) % 3;
                                GetNeighbor(neighbor, orient, out neighbor2, out nid2);
                            }

                            Indices[neighbor.ID, orient] = count;
                        }

                        count++;
                    }
                }
            }
        }

        private void GetNeighbor(ITriangle tri, int i, out ITriangle neighbor, out int nid)
        {
            neighbor = tri.GetNeighbor(i);
            nid = neighbor == null ? -1 : neighbor.ID;
        }
    }
}