Changed type of Mesh.Points from vector<Point*> to vector<Point>

Most methods in Edge and Mesh were altered.
Methods in Edge which required access to Mesh.Points were moved into Mesh

src/Mesh/src/Edge.cpp

 #include "Mesh.hpp"
 
-Edge::Edge(Curve *c, bool direction) {
+Edge::Edge(Curve *c, bool direction, size_t v) {
     ConstraintCurve = c;
-    Node = new Point(direction ? c->start() : c->end());
+    Node = v;
     Twin = this;
     Next = nullptr;
     Prev = nullptr;
     Orientation = direction;
 }
 
-Edge::Edge(Curve *c, bool direction, Point const *v) {
-    ConstraintCurve = c;
-    Node = v; //Should have C->point(0.0) == *(v) || c->point(1.0) == *(v)
-    Twin = this;
-    Next = nullptr;
-    Prev = nullptr;
-    Orientation = direction;
-}
-
-bool Edge::is_protruding() const {
-    /*
-        See Chapter 1.4 of "Triangulations and Applications" by Øyvind Hjelle and Morten Dæhlen
-    */
-    Point const *v0 = Prev->Node;
-    Point const *v1 = Node;
-    Point const *v2 = Next->Node;
-
-    double v1x = v2->X - v1->X;
-    double v1y = v2->Y - v1->Y;
-
-    double v0x = v0->X - v1->X;
-    double v0y = v0->Y - v1->Y;
-
-    double area = v1x * v0y - v1y * v0x;
-
-    if (area > 0.0) {
-        // Make sure no boundary points interior to triangle
-        // Calculate barrycentric coordinates of p
-        Edge * e = Next->Next;
-
-        while (e != Prev) {
-            Point const *p = e->Node;
-
-            double v2x = v2->X - p->X;
-            double v2y = v2->Y - p->Y;
-
-            double v1x = v1->X - p->X;
-            double v1y = v1->Y - p->Y;
-
-            double v0x = v0->X - p->X;
-            double v0y = v0->Y - p->Y;
-
-            double b0 = (v0x * v1y - v0y * v1x);
-            double b1 = (v1x * v2y - v1y * v2x);
-            double b2 = (v2x * v0y - v2y * v0x);
-
-            if (b0 >= 0.0 && b0 <= area && b1 >= 0.0 && b1 <= area && b2 >= 0.0 && b2 <= area) {
-                return false;    // Point is interior to triangle
-            } else {
-                e = e->Next;
-            }
-        }
-        return true;
-    } else {
-        return false;
-    }
-}
-
-bool Edge::is_locally_optimal() const {
-    /*
-        See Chapter 3.7 of "Triangulations and Applications" by Øyvind Hjelle and Morten Dæhlen
-    */
-
-    if (ConstraintCurve != nullptr) {
-        return true;
-    } else {
-        Point const *p3 = Node;
-        Point const *p2 = Prev->Node;
-        Point const *p1 = Twin->Node;
-        Point const *p4 = Twin->Prev->Node;
-
-        double v1x = p3->X - p2->X;
-        double v1y = p3->Y - p2->Y;
-        double v2x = p1->X - p2->X;
-        double v2y = p1->Y - p2->Y;
-        double v3x = p1->X - p4->X;
-        double v3y = p1->Y - p4->Y;
-        double v4x = p3->X - p4->X;
-        double v4y = p3->Y - p4->Y;
-
-        double d1 = sqrt(v1x * v1x + v1y * v1y);
-        double d2 = sqrt(v2x * v2x + v2y * v2y);
-        double d3 = sqrt(v3x * v3x + v3y * v3y);
-        double d4 = sqrt(v4x * v4x + v4y * v4y);
-        double tol = -(d1 * d2 * d3 * d4 * FLT_EPSILON);
-
-        double sina = v1x * v2y - v1y * v2x;
-        double sinb = v3x * v4y - v3y * v4x;
-        double cosa = v1x * v2x + v1y * v2y;
-        double cosb = v3x * v4x + v3y * v4y;
-        double cct = sina * cosb + cosa * sinb;
-
-        return (cct < tol ? false : true);
-    }
-}
-
 bool Edge::is_valid() const {
     bool value = true;
 
@@ -115,70 +19,6 @@ bool Edge::is_valid() const {
     return value;
 }
 
-bool Edge::is_encroached(Point const *v2) const {
-    /*
-    A constrained edge is encroached if a triangle and it's circumcenter lie on opposite sides of the edge.
-    This is equivalent to a node being in the diameteral ball of the edge?
-    This only occurs if one of the triangles attached to the edge encroaches the edge?
-    This only occurs if the angle of the triangles attached to the edge has an angle opposite the edge of greater than 90 degrees.
-    Using the dot product, this requires that the normalized dot product < cos(90) = 0
-    */
-
-    if (ConstraintCurve == nullptr) {
-        return false;
-    } else {
-        Point const *v0 = base();
-        Point const *v1 = tip();
-
-        double dx0 = v0->X - v2->X;
-        double dy0 = v0->Y - v2->Y;
-
-        double dx1 = v1->X - v2->X;
-        double dy1 = v1->Y - v2->Y;
-
-        double dot = dx0 * dx1 + dy0 * dy1;
-        double tol = sqrt(dx0 * dx0 + dy0 * dy0) * sqrt(dx1 * dx1 + dy1 * dy1) * FLT_EPSILON;
-
-        return (dot < tol ? true : false);
-    }
-}
-
-bool Edge::is_attached(Point const *p, Edge *&e) const {
-    if (*this->tip() == *p) {
-        return true;
-    }
-
-    e = const_cast<Edge *>(this);
-
-    if (e != e->Twin) {
-        e = e->Twin->Next;
-        while (e != this) {
-            if (*e->tip() == *p) {
-                return true;
-            } else if (e->Twin != e) {
-                e = e->Twin->Next;
-            } else {
-                break;
-            }
-        }
-    }
-
-    if (e == e->Twin) {
-        e = this->Prev;
-        while (e != e->Twin) {
-            if (*e->base() == *p) {
-                e = e->Twin;
-                return true;
-            } else {
-                e = e->Twin->Prev;
-            }
-        }
-    }
-
-    return false;
-}
-
-// Topology Altering Methods
 bool Edge::swap() {
     if (ConstraintCurve == nullptr) {
         Edge * e1 = Next;
@@ -219,141 +59,6 @@ bool Edge::swap() {
     }
 }
 
-bool Edge::recursive_swap() {
-    // TODO, May need to have two different recursive swap methods, one for midpoint insertion and one for circumcenter insertion
-    if (!is_locally_optimal() && swap()) {
-        Edge * enext = Next;
-        Edge * eprev = Prev;
-        Edge * tnext = Twin->Next;
-        Edge * tprev = Twin->Prev;
-
-        enext->recursive_swap();
-        eprev->recursive_swap();
-        tnext->recursive_swap();
-        tprev->recursive_swap();
-
-        return true;
-    } else {
-        return false;
-    }
-}
-
-void Edge::split_edge(std::vector<Point const *> &points, std::vector<Edge *> &edges) {
-    /*
-        Splits edge into two edges at the midpoint without creating any new triangles.
-        Used for initial polygon refinement.
-    */
-
-    // TODO: Refactor into split_constrainted_edge
-
-    Point * p = new Point;
-    points.push_back(p);
-
-    Curve *c;
-    if (ConstraintCurve != nullptr) { // Constrained Edge
-        Vertex *v = new Vertex;    // TODO: Need to manage this memory. Destructors.
-        c = ConstraintCurve->split(v, 0.5);
-
-        p->X = v->x();
-        p->Y = v->y();
-    } else { // Unconstrained Edge
-        c = nullptr;
-        *(p) = Point((Node->X + Next->Node->X) / 2.0, (Node->Y + Next->Node->Y) / 2.0);
-    }
-
-    if (this == Twin) { // Boundary Edge
-        Edge * e = new Edge;
-        edges.push_back(e);
-
-        // Constraint Curve
-        e->Orientation = Orientation;
-        if (Orientation) {
-            e->ConstraintCurve = c;
-        } else {
-            e->ConstraintCurve = ConstraintCurve;
-            ConstraintCurve = c;
-        }
-
-        // Connectivity
-        e->Node = p;
-        e->Next = Next;
-        e->Prev = this;
-        e->Twin = e;
-        e->Mark = false;
-
-        Next->Prev = e;
-        Next->Mark = false;
-
-        Next = e;
-        Mark = false;
-    } else {
-        throw std::exception(); // Function should only be called on
-    }
-}
-
-Point Edge::circumcenter() const {
-    double xa = Node->X;
-    double ya = Node->Y;
-    double xb = Prev->Node->X - xa;
-    double yb = Prev->Node->Y - ya;
-    double xc = Next->Node->X - xa;
-    double yc = Next->Node->Y - ya;
-
-    double d = 2.0 * (xb * yc - yb * xc);
-    double db = xb * xb + yb * yb;
-    double dc = xc * xc + yc * yc;
-    double x = (yc * db - yb * dc) / d + xa;
-    double y = (xb * dc - xc * db) / d + ya;
-
-    return Point{x, y};
-}
-
-double Edge::circumradius() const {
-    double xa = Node->X;
-    double ya = Node->Y;
-    double xb = Prev->Node->X - xa;
-    double yb = Prev->Node->Y - ya;
-    double xc = Next->Node->X - xa;
-    double yc = Next->Node->Y - ya;
-    xa = xb - xc;
-    ya = yb - yc;
-
-    double den = 2.0 * abs(xb * yc - yb * xc);
-    double num = sqrt(xa * xa + ya * ya) * sqrt(xb * xb + yb * yb) * sqrt(xc * xc + yc * yc);
-
-    return num / den;
-}
-
-double Edge::length() const {
-    double dx = Node->X - Next->Node->X;
-    double dy = Node->Y - Next->Node->Y;
-
-    return sqrt(dx * dx + dy * dy);
-}
-
-double Edge::shortest_edge_length() const {
-    double x0 = Node->X;
-    double y0 = Node->Y;
-    double x1 = Next->Node->X;
-    double y1 = Next->Node->Y;
-    double x2 = Prev->Node->X;
-    double y2 = Prev->Node->Y;
-
-    double dx = x0 - x1;
-    double dy = y0 - y1;
-    double dl = dx * dx + dy * dy;
-
-    dx = x1 - x2;
-    dy = y1 - y2;
-    dl = fmin(dl, dx * dx + dy * dy);
-
-    dx = x2 - x0;
-    dy = y2 - y0;
-    dl = fmin(dl, dx * dx + dy * dy);
-
-    return sqrt(dl);
-}
-
 void Edge::recursive_mark() {
     if (Mark && Next->Mark && Prev->Mark) {
         Next->Mark = false;

src/Mesh/src/Edge.h

 #ifndef OERSTED_EDGE_H
 #define OERSTED_EDGE_H
 
-#include "Mesh.h"
+#include "Point.h"
+
+class Mesh;
 
 class Edge {
 public:
     friend class Mesh;
 
-    friend bool are_intersecting(Edge const *e0, Edge const *e1);
+    friend bool are_intersecting(Edge const *e0, Edge const *e1, Mesh const &m);
 
     // Constructors
-    Edge() : Node(nullptr), Next(nullptr), Prev(nullptr), Twin(this), ConstraintCurve(nullptr), Orientation(true), Mark(false) {};
-
-    Edge(Point &v, Edge &n, Edge &p, Edge &t) : Node(&v), Next(&n), Prev(&p), Twin(&t), ConstraintCurve(nullptr), Orientation(true), Mark(false) {};
+    Edge() : Node(SIZE_MAX), Next(nullptr), Prev(nullptr), Twin(this), ConstraintCurve(nullptr), Orientation(true), Mark(false) {};
 
-    Edge(Curve *c, bool Orientation);
+    Edge(size_t v, Edge &n, Edge &p, Edge &t) : Node(v), Next(&n), Prev(&p), Twin(&t), ConstraintCurve(nullptr), Orientation(true), Mark(false) {};
 
-    Edge(Curve *c, bool Orientation, Point const *v);
+    Edge(Curve *c, bool Orientation, size_t v);
 
     // Accessors
-    Point const *node() const { return Node; };
+    size_t node() const { return Node; };
 
-    Point const *base() const { return Node; };
+    size_t base() const { return Node; };
 
-    Point const *tip() const { return (next() == nullptr ? twin()->base() : next()->base()); };
+    size_t tip() const { return (next() == nullptr ? twin()->base() : next()->base()); };
 
     Edge const *next() const { return Next; };
 
@@ -41,36 +41,16 @@ public:
         return (node() == e.node()) && (constraint_curve() == e.constraint_curve()) && (twin()->node() == e.twin()->node());
     };
 
-    Point circumcenter() const;
-
-    double circumradius() const;
-
-    double length() const;
-
-    double shortest_edge_length() const;
-
-    bool is_protruding() const;
-
     bool is_constrained() const { return (ConstraintCurve != nullptr); };
 
-    bool is_locally_optimal() const;
-
     bool is_valid() const;
 
-    bool is_encroached(Point const *p) const;
-
-    bool is_attached(Point const *p, Edge *&e) const;
-
     bool swap();
 
-    bool recursive_swap();
-
-    void split_edge(std::vector<Point const *> &verts, std::vector<Edge *> &edges);
-
     void recursive_mark();
 
 protected:
-    Point const *Node;      //Start of edge
+    size_t Node;            //Start of edge
     Edge *Next;             //In triangle
     Edge *Twin;             //Adjacent triangle
     Edge *Prev;             //In triangle

src/Mesh/src/Mesh.h

@@ -2,6 +2,8 @@
 #define OERSTED_MESH_H
 
 #include "Sketch.hpp"
+#include "Point.h"
+#include "Edge.h"
 
 #include <fstream>
 #include <algorithm>
@@ -15,11 +17,9 @@ enum class InsertPointResult {
     Success, Midpoint, Duplicate, Failure
 };
 
-class Point;
-
-class Edge;
-
 class Mesh {
+    friend class Edge;
+
 public:
     double MinimumElementQuality = 0.0;
     double MinimumElementSize = 0.0;
@@ -31,7 +31,9 @@ public:
 
     void create();
 
-    Point const *point(size_t i) const { return Points[i]; };
+    Point circumcenter(Edge const *e) const;
+
+    Point const point(size_t i) const { return Points[i]; };
 
     Edge const *edge(size_t i) const { return Edges[i]; };
 
@@ -51,20 +53,18 @@ public:
 
     void save_as(std::string path, std::string file_name) const;
 
-    LocateTriangleResult locate_triangle(Point const *p, Edge *&e) const;
+    LocateTriangleResult locate_triangle(Point const &p, Edge *&e) const;
 
-    LocateTriangleResult locate_triangle(Point const *p, Edge const *&e) const;
+    LocateTriangleResult locate_triangle(Point const &p, Edge const *&e) const;
 
-    LocateTriangleResult locate_triangle(Point const *p) const {
+    LocateTriangleResult locate_triangle(Point const p) const {
         Edge const *e = Edges.back();
         return locate_triangle(p, e);
     };
 
-    bool in_triangle(Point const *p, Edge const *&e) const;
+    InsertPointResult insert_point(Point const &p, Edge *e);
 
-    InsertPointResult insert_point(Point const *p, Edge *e);
-
-    InsertPointResult insert_point(Point const *p) { return insert_point(p, Edges.back()); };
+    InsertPointResult insert_point(Point const p) { return insert_point(p, Edges.back()); };
 
     InsertPointResult insert_circumcenter(Edge *e);
 
@@ -76,13 +76,31 @@ public:
 
     void refine_once(std::vector<size_t> index, std::vector<double> circumradius, std::vector<double> quality);
 
-    bool edges_are_valid();
+    bool edges_are_valid() const;
+
+    double circumradius(Edge const *e) const;
+
+    double length(Edge const *e) const;
+
+    double shortest_edge_length(Edge const *e) const;
+
+    bool is_protruding(Edge const *e) const;
+
+    bool recursive_swap(Edge *edge) const;
+
+    bool is_locally_optimal(Edge const *edge) const;
+
+    void split_edge(Edge *edge);
+
+    bool is_encroached(Edge const *e, Point const p2) const;
+
+    bool is_attached(Edge *&edge_out, Point const &p) const;
 
 protected:
     Contour const *Boundary;
     std::vector<Curve const *> Curves;
     std::vector<Contour const *> Contours;
-    std::vector<Point const *> Points;
+    std::vector<Point> Points;
     std::vector<Edge *> Edges;
     std::vector<Edge *> Triangles;
 

src/Mesh/src/util.cpp

 #include "Mesh.hpp"
 
-bool are_intersecting(const Edge *e0, const Edge *e1) {
+// TODO: Make these mesh routines
+bool are_intersecting(Edge const *e0, Edge const *e1, Mesh const &m) {
     // TODO, Make more detailed return type enumeration
     if (e0->ConstraintCurve != nullptr && e0->ConstraintCurve == e1->ConstraintCurve) {
         return false;
     }
 
-    const Point *v00 = e0->base();
-    const Point *v01 = e0->tip();
-    const Point *v10 = e1->base();
-    const Point *v11 = e1->tip();
+    Point const v00 = m.point(e0->base());
+    Point const v01 = m.point(e0->tip());
+    Point const v10 = m.point(e1->base());
+    Point const v11 = m.point(e1->tip());
 
-    double xs0 = (v00->X + v01->X) / 2.0;
-    double ys0 = (v00->Y + v01->Y) / 2.0;
-    double xs1 = (v10->X + v11->X) / 2.0;
-    double ys1 = (v10->Y + v11->Y) / 2.0;
+    double xs0 = (v00.X + v01.X) / 2.0;
+    double ys0 = (v00.Y + v01.Y) / 2.0;
+    double xs1 = (v10.X + v11.X) / 2.0;
+    double ys1 = (v10.Y + v11.Y) / 2.0;
 
-    double xd0 = (v00->X - v01->X) / 2.0;
-    double yd0 = (v00->Y - v01->Y) / 2.0;
-    double xd1 = (v10->X - v11->X) / 2.0;
-    double yd1 = (v10->Y - v11->Y) / 2.0;
+    double xd0 = (v00.X - v01.X) / 2.0;
+    double yd0 = (v00.Y - v01.Y) / 2.0;
+    double xd1 = (v10.X - v11.X) / 2.0;
+    double yd1 = (v10.Y - v11.Y) / 2.0;
 
     double d0 = xd0 * xd0 + yd0 * yd0;
     double d1 = xd1 * xd1 + yd1 * yd1;
@@ -27,10 +28,9 @@ bool are_intersecting(const Edge *e0, const Edge *e1) {
     double tol = (d0 * d1) * FLT_EPSILON;
 
     if (cross < tol) {
-        /*
-            Lines are nearly parallel
-            There are four possible minimum distance points between the lines
-        */
+        // Lines are nearly parallel
+        // There are four possible minimum distance points between the lines
+
         double s, dx, dy, dmin = DBL_MAX;
 
         s = ((xd0 - xd1) * (xs0 - xs1) + (yd0 - yd1) * (ys0 - ys1)) /
@@ -59,8 +59,7 @@ bool are_intersecting(const Edge *e0, const Edge *e1) {
 
         tol = (d0 + d1) * FLT_EPSILON;
         return (dmin < tol);
-    } else {
-        // Lines are not parallel
+    } else { // Lines are not parallel
         double s0 = abs(xd1 * (ys0 - ys1) - yd1 * (xs0 - xs1));
         double s1 = abs(xd0 * (ys0 - ys1) - yd0 * (xs0 - xs1));
         tol = cross * (1.0 - FLT_EPSILON);
@@ -69,23 +68,66 @@ bool are_intersecting(const Edge *e0, const Edge *e1) {
     }
 }
 
-void element_quality(std::vector<Edge *> &triangle, std::vector<double> &radii, std::vector<double> &quality) {
+bool in_triangle(Point const &p, Edge const *&e, Mesh const &m) {
+    double xp = p.X;
+    double yp = p.Y;
+
+    Point const p0 = m.point(e->node());
+    Point const p1 = m.point(e->next()->node());
+    Point const p2 = m.point(e->prev()->node());
+
+    double dx0p = p0.X - xp;
+    double dy0p = p0.Y - yp;
+
+    double dx1p = p1.X - xp;
+    double dy1p = p1.Y - yp;
+
+    double dx2p = p2.X - xp;
+    double dy2p = p2.Y - yp;
+
+    double dx01 = p0.X - p1.X;
+    double dy01 = p0.Y - p1.Y;
+
+    double dx12 = p1.X - p2.X;
+    double dy12 = p1.Y - p2.Y;
+
+    double dx20 = p2.X - p0.X;
+    double dy20 = p2.Y - p0.Y;
+
+    double area012 = dx01 * dy12 - dy01 * dx12;
+
+    double tol = FLT_EPSILON * area012;
+
+    double area01p = dx0p * dy1p - dx1p * dy0p;
+    double area12p = dx1p * dy2p - dx2p * dy1p;
+    double area20p = dx2p * dy0p - dx0p * dy2p;
+
+    return (area01p > -tol && area12p > -tol && area20p > -tol);
+}
+
+void element_quality(std::vector<Edge *> &triangle, std::vector<double> &radii, std::vector<double> &quality, Mesh const &m) {
     radii.resize(0);
     quality.resize(0);