Mesh construction now determines which contour contains each triangle

Paves the way for local region refinement. Necessary for converting refineable mesh to FEA discretization.

src/Mesh/CMakeLists.txt

@@ -9,6 +9,7 @@ set(SOURCE_FILES
         ./src/InsertionQueuer.h ./src/InsertionQueuer.cpp
         ./src/Mesh.h    ./src/Mesh.cpp
         ./src/Point.h   ./src/Point.cpp
+        ./src/Triangle.h ./src/Triangle.cpp
         )
 
 add_library(mesh SHARED ${SOURCE_FILES})

src/Mesh/src/Edge.h

@@ -31,6 +31,8 @@ public:
 
     size_t constraint() const { return Constraint; };
 
+    bool is_constrained() const {return Constraint != 0; };
+
     bool orientation() const { return Orientation; };
 
     bool mark() const { return Mark; };

src/Mesh/src/Mesh.cpp

@@ -541,7 +541,6 @@ size_t Mesh::num_edges() const {
 void Mesh::create() {
     create_boundary_polygon();
     triangulate_boundary_polygon();
-    get_triangles();
 
     insert_internal_boundaries();
     get_triangles();
@@ -603,8 +602,8 @@ void Mesh::element_quality(std::vector<double> &radii, std::vector<double> &qual
     radii.reserve(Triangles.size());
     quality.reserve(Triangles.size());
     for (size_t i = 0; i != Triangles.size(); ++i) {
-        double r = circumradius(Triangles[i]);
-        double l = shortest_edge_length(Triangles[i]);
+        double r = circumradius(Triangles[i].Edge);
+        double l = shortest_edge_length(Triangles[i].Edge);
 
         radii.push_back(r);
         quality.push_back(l / r / sqrt(3.0)); // sqrt(3.0) = (shortest edges length) / radius of equilateral triangle
@@ -615,10 +614,43 @@ void Mesh::get_triangles() {
     Triangles.resize(0);
     Triangles.reserve(2 * num_points());
 
-    mark_triangles();
     for (auto &e : Edges) {
-        if (e.Mark) {
-            Triangles.push_back(e.Self);
+        e.Mark = true;
+    }
+
+    std::vector<size_t> edge_queue;
+    std::vector<size_t> contours;
+    edge_queue.reserve(2 * num_points());
+    contours.reserve(2 * num_points());
+    for (size_t contour_index = 0; contour_index != Contours.size(); ++contour_index) {
+        edge_queue.resize(0);
+
+        std::shared_ptr<Contour const> &contour = Contours[contour_index];
+        std::shared_ptr<Curve const> c = contour->curve(0);
+        std::shared_ptr<BoundaryConstraint> bc = boundary_constraint(c);
+        DartConstraint const dc = dart_constraint(bc->dart(0));
+
+        edge_queue.push_back(contour->orientation(0) ? dc.forward_dart() : dc.reverse_dart());
+        while(edge_queue.size() > 0) {
+            Edge &et = Edges[edge_queue.back()];
+            Edge &en = Edges[et.next()];
+            Edge &ep = Edges[et.prev()];
+
+            edge_queue.pop_back();
+
+            if (et.Mark) {
+                Triangles.emplace_back(et.self(), contour_index);
+
+                if(!en.is_constrained() && en.Mark) {
+                    edge_queue.push_back(en.twin());
+                    en.Mark = false;
+                }
+
+                if(!ep.is_constrained() && ep.Mark) {
+                    edge_queue.push_back(ep.twin());
+                    ep.Mark = false;
+                }
+            }
         }
     }
 }
@@ -765,7 +797,7 @@ void Mesh::refine_once(std::vector<size_t> index, std::vector<double> radii, std
     for (size_t i = 0; i != Triangles.size(); ++i) {
         size_t j = index[i];
         if ((triangle(j).Mark) && ((radii[j] > MaximumElementSize) || (radii[j] > MinimumElementSize && quality[j] < MinimumElementQuality))) {
-            add_circumcenter_to_queue(Triangles[j]);
+            add_circumcenter_to_queue(Triangles[j].Edge);
             insert_from_queue();
         }
     }

src/Mesh/src/Mesh.h

@@ -12,6 +12,7 @@
 #include "Edge.h"
 #include "InsertionQueuer.h"
 #include "Point.h"
+#include "Triangle.h"
 
 enum class LocateTriangleResult {
     Interior, Exterior, Boundary, Point // TODO: Enumerate cases in function when triangle is located by a point near the boundary
@@ -21,7 +22,7 @@ enum class InsertPointResult {
     Success, Midpoint, Duplicate, Failure
 };
 
-class Mesh {
+class Mesh { // TODO: Namespaces. Also, need to rename a bunch of things to be more consistent and clear w.r.t. how various arrays are accessed (e.g. by edge/dart index versus plain index)
 public:
     Mesh(Sketch &s);
 
@@ -131,7 +132,7 @@ public:
 
     Edge const &twin(Edge const &e) const { return Edges[e.Twin]; };
 
-    Edge const &triangle(size_t i) const { return Edges[Triangles[i]]; };
+    Edge const &triangle(size_t i) const { return Edges[Triangles[i].Edge]; }; // TODO: Rename to triangle_dart
 
     LocateTriangleResult locate_triangle(Point const p, size_t &ei) const;
 
@@ -158,7 +159,7 @@ protected:
 
     std::vector<std::shared_ptr<BoundaryConstraint>> BoundaryConstraints;
 
-    std::vector<size_t> Triangles;
+    std::vector<Triangle> Triangles;
 
     std::vector<std::unique_ptr<InsertionQueuer>> Queue;
 

src/Mesh/src/Triangle.cpp

+#include "Triangle.h"

src/Mesh/src/Triangle.h

+#ifndef OERSTED_TRIANGLE_H
+#define OERSTED_TRIANGLE_H
+
+#include <cstddef>
+#include <cstdint>
+
+class Triangle {
+public:
+    Triangle() : Edge{SIZE_MAX}, Contour{SIZE_MAX} {};
+    Triangle(size_t e, size_t c) : Edge{e}, Contour{c} {};
+
+    size_t Edge;
+    size_t Contour;
+};
+
+#endif //OERSTED_TRIANGLE_H