Fix for inserting circumcenter of triangles which lie outside the discretized domain

Modification was necessary to introduce uniform boundary discretizations.

More work may be needed to make this fully robust.

src/Mesh/src/BoundaryConstraint.cpp

@@ -33,6 +33,33 @@ void BoundaryConstraint::refine(Mesh &m, double_t tol) {
             }
         }
 
+        m.insert_from_queue();
+    }
+}
+
+void BoundaryConstraint::make_uniform(Mesh &m) {
+    double_t delta_min{DBL_MAX};
+
+    for (size_t i : DartConstraints) {
+        DartConstraint const &dc = m.dart_constraint(i);
+        delta_min = std::min(delta_min, std::abs(dc.S1 - dc.S0));
+    }
+    delta_min *= (1 + FLT_EPSILON);
+
+    double_t delta_max{DBL_MAX};
+    while (delta_max > delta_min) {
+        delta_max = delta_min;
+
+        for(size_t i : DartConstraints) {
+            DartConstraint const &dc = m.dart_constraint(i);
+            double_t delta = std::abs(dc.S1 - dc.S0);
+            if (delta > delta_min) {
+                m.add_to_queue(std::make_unique<MidpointQueuer>(dc.dart()));
+
+                delta_max = std::max(delta_max, delta);
+            }
+        }
+
         m.insert_from_queue();
     }
 }
\ No newline at end of file

src/Mesh/src/BoundaryConstraint.h

@@ -26,9 +26,11 @@ public:
 
     void refine(Mesh &m, double_t tol);
 
-    size_t size_dart_constraints() const {return DartConstraints.size(); };
+    void make_uniform(Mesh &m);
 
-    size_t dart(size_t i) const {return DartConstraints[i]; };
+    size_t size_dart_constraints() const { return DartConstraints.size(); };
+
+    size_t dart(size_t i) const { return DartConstraints[i]; };
 
 protected:
     std::vector<size_t> DartConstraints;

src/Mesh/src/Mesh.cpp

@@ -748,6 +748,12 @@ void Mesh::insert_internal_boundaries() {
         */
     }
 
+    /*
+    for (auto bc : BoundaryConstraints) {
+        bc->make_uniform(*this);
+    }
+    */
+
     split_encroached_edges();
 
     // TODO: Enforce BoundaryConstraint conditions (e.g. uniform discretization, boundary maps)
@@ -944,6 +950,7 @@ void Mesh::triangulate_boundary_polygon() {
 }
 
 LocateTriangleResult Mesh::locate_triangle(Point const p, size_t &ei) const {
+    // TODO: Check for visibility of vertex? That someone contradicts the purpose of this method as a general point location function (as opposed to purely for mesh refinement)
     double xp = p.X;
     double yp = p.Y;
 
@@ -1314,13 +1321,20 @@ InsertPointResult Mesh::insert_midpoint(size_t ei) {
     return InsertPointResult::Midpoint;
 }
 
-AddToQueueResult Mesh::add_point_to_queue(Point const vc, size_t ei) {
+AddToQueueResult Mesh::add_point_to_queue(Point vc, size_t ei) {
     // Find triangle containing point
     LocateTriangleResult result = locate_triangle(vc, ei);
 
-    // Do not insert point if it is duplicate
     if (result == LocateTriangleResult::Point) {
-        return AddToQueueResult::Duplicate;
+        return AddToQueueResult::Duplicate; // Do not insert point if it is duplicate
+    }
+
+    if (result == LocateTriangleResult::Exterior) {
+        if (is_constrained(ei)) { // Move attempted insertion point to the midpoint of the located edge to force encroachement
+            vc = midpoint(ei);
+        } else {
+            throw std::exception(); // TODO: No triangle found containing circumcenter, boundary edge was encroached by circumcenter, and the located edge was not a constrained edge
+        }
     }
 
     // Test edges in current and adjacent triangles for encroachment

src/Mesh/src/Mesh.h

@@ -161,27 +161,36 @@ public:
 
     DartConstraint const dart_constraint(size_t i) { return DartConstraints[i]; };
 
-    Point const base(Edge const e) const { return Points[e.Node]; };
+    Point const &base(Edge const &e) const { return Points[e.Node]; };
 
-    Point const base(size_t ei) const { return Points[node(ei)]; };
+    Point const &base(size_t i) const { return Points[node(i)]; };
 
-    Point const point(size_t i) const { return Points[i]; };
+    Point const &point(size_t i) const { return Points[i]; };
 
-    Point const point(Edge const e) const { return Points[e.Node]; };
+    Point const &point(Edge const &e) const { return Points[e.Node]; };
 
-    Point const tip(Edge const e) const { return Points[next(e).Node]; };
+    Point const &tip(Edge const &e) const { return Points[next(e).Node]; };
 
-    Point const tip(size_t ei) const { return Points[node(next(ei))]; };
+    Point const &tip(size_t i) const { return Points[node(next(i))]; };
 
-    Edge const edge(size_t i) const { return Edges[i]; };
+    Point midpoint(size_t e) { return midpoint(Edges[e]); };
 
-    Edge const next(Edge const e) const { return Edges[e.Next]; };
+    Point midpoint(Edge const e) const {
+        Point const &p0 = base(e);
+        Point const &p1 = tip(e);
 
-    Edge const prev(Edge const e) const { return Edges[e.Prev]; };
+        return Point((p0.X + p1.X) / 2.0, (p0.Y + p1.Y) / 2.0);
+    }
+
+    Edge const &edge(size_t i) const { return Edges[i]; };
+
+    Edge const &next(Edge const &e) const { return Edges[e.Next]; };
+
+    Edge const &prev(Edge const &e) const { return Edges[e.Prev]; };
 
-    Edge const twin(Edge const e) const { return Edges[e.Twin]; };
+    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]]; };
 
     LocateTriangleResult locate_triangle(Point const p, size_t &ei) const;
 

test/UseCases/test_Rotor.cpp

@@ -46,8 +46,8 @@ TEST(Rotor, Suite0) {
 
     Mesh m{s};
     m.MaximumElementSize = 2.0;
-    m.MinimumElementSize = 0.20;
-    m.MinimumElementQuality = 0.1;
+    m.MinimumElementSize = 0.02;
+    m.MinimumElementQuality = 0.3;
 
     m.create();
 
@@ -140,8 +140,8 @@ TEST(Rotor, Circular_Barrier_Syncrel) {
     // Mesh
     Mesh m{s};
     m.MaximumElementSize = 1.0;
-    m.MinimumElementSize = 0.1;
-    m.MinimumElementQuality = 0.1;
+    m.MinimumElementSize = 0.01;
+    m.MinimumElementQuality = 0.3;
 
     m.create();
 

test/UseCases/test_Stator.cpp

@@ -73,9 +73,9 @@ TEST(Stator, 0) {
 
     Mesh mesh{sketch};
 
-    mesh.MinimumElementQuality = 0.1;
+    mesh.MinimumElementQuality = 0.3;
     mesh.MaximumElementSize = 4.0;
-    mesh.MinimumElementSize = 0.40;
+    mesh.MinimumElementSize = 0.040;
 
     mesh.create();
 
@@ -228,9 +228,9 @@ TEST(Stator, 1) {
     // Create Mesh
     Mesh mesh{sketch};
 
-    mesh.MinimumElementQuality = 0.1;
+    mesh.MinimumElementQuality = 0.3;
     mesh.MaximumElementSize = M_PI * radius_1_5_0->dim() / Nt / 3.0;
-    mesh.MinimumElementSize = mesh.MaximumElementSize / 10.0;
+    mesh.MinimumElementSize = mesh.MaximumElementSize / 100.0;
 
     mesh.create();