Add length based boundary curve refinement to initial mesh refinement

src/Mesh/src/BoundaryConstraint.cpp

 #include "BoundaryConstraint.h"
 #include "DartConstraint.h"
+#include "Mesh.h"
 
 void BoundaryConstraint::add_to_queue(std::vector<std::unique_ptr<InsertionQueuer>> &queue, std::vector<DartConstraint> &constraints, size_t dart_constraint) {
     if (UniformDiscretization) {
@@ -11,4 +12,27 @@ void BoundaryConstraint::add_to_queue(std::vector<std::unique_ptr<InsertionQueue
         size_t dart = constraints[dart_constraint].dart();
         queue.push_back(std::make_unique<MidpointQueuer>(dart));
     }
+}
+
+void BoundaryConstraint::refine(Mesh &m, double_t tol) {
+    double_t err_max{DBL_MAX};
+
+    while (err_max > tol) {
+        err_max = 0.0;
+        for (size_t i : DartConstraints) {
+            DartConstraint const &dc = m.dart_constraint(i);
+
+            double_t len_c = ConstraintCurve->length(dc.S0, dc.S1);
+            double_t len_e = m.edge(dc.dart()).length(m);
+            double_t err_i = std::abs(len_c - len_e) / len_c;
+
+            if (err_i > tol) {
+                m.add_to_queue(std::make_unique<MidpointQueuer>(dc.dart()));
+
+                err_max = std::max(err_max, err_i);
+            }
+        }
+
+        m.insert_from_queue();
+    }
 }
\ No newline at end of file

src/Mesh/src/BoundaryConstraint.h

@@ -24,6 +24,8 @@ public:
 
     void add_to_queue(std::vector<std::unique_ptr<InsertionQueuer>> &queue, std::vector<DartConstraint> &constraints, size_t dart);
 
+    void refine(Mesh &m, double_t tol);
+
     size_t size_dart_constraints() const {return DartConstraints.size(); };
 
     size_t dart(size_t i) const {return DartConstraints[i]; };

src/Mesh/src/Edge.cpp

-#include "Mesh.hpp"
\ No newline at end of file
+#include "Edge.h"
+#include "Mesh.h"
+
+double_t Edge::length(Mesh const &m) const {
+    return dist(m.point(node()), m.point(m.edge(next()).node()));
+}
\ No newline at end of file

src/Mesh/src/Edge.h

@@ -19,10 +19,6 @@ public:
 
     Edge(size_t n, size_t s, size_t c, bool d) : Node(n), Self(s), Next(SIZE_MAX), Twin(SIZE_MAX), Prev(SIZE_MAX), Constraint(c), Orientation(d), Mark(false) {};
 
-    void add_to_queue(std::vector<std::unique_ptr<InsertionQueuer>> &queue, std::vector<DartConstraint> &constraints) {
-        constraints[Constraint].add_to_queue(queue, constraints);
-    }
-
     size_t node() const { return Node; };
 
     size_t self() const { return Self; };
@@ -49,6 +45,12 @@ public:
                 Orientation == e.Orientation);
     };
 
+    double_t length(Mesh const &m) const;
+
+    void add_to_queue(std::vector<std::unique_ptr<InsertionQueuer>> &queue, std::vector<DartConstraint> &constraints) {
+        constraints[Constraint].add_to_queue(queue, constraints);
+    }
+
 protected:
     size_t Node;            //Point at start of this edge
     size_t Self;            //This edge in this triangle

src/Mesh/src/Mesh.cpp

@@ -738,6 +738,7 @@ void Mesh::insert_internal_boundaries() {
     // TODO: Do length based refinement instead
     //      There are some errors that occur with half-circles becoming degenerate when the LineSegment is not encroached because this isn't done
     for (auto bc : BoundaryConstraints) {
+        bc->refine(*this, 0.01);
         /*
         auto &temp = *(bc.get());
         if(bc->size_dart_constraints() == 1) {

src/Mesh/src/Mesh.h

@@ -8,9 +8,10 @@
 
 #include "Sketch.hpp"
 
-#include "Point.h"
+#include "BoundaryConstraint.h"
 #include "Edge.h"
 #include "InsertionQueuer.h"
+#include "Point.h"
 
 enum class LocateTriangleResult {
     Interior, Exterior, Boundary, Point // TODO: Enumerate cases in function when triangle is located by a point near the boundary
@@ -86,8 +87,8 @@ class Mesh {
 public:
     Mesh(Sketch &s);
 
+    friend class BoundaryConstraint;
     friend class CircumcenterQueuer;
-
     friend class MidpointQueuer;
 
     double MinimumElementQuality = 0.0;
@@ -146,6 +147,8 @@ public:
 
     size_t twin(size_t ei) const { return Edges[ei].Twin; };
 
+    void add_to_queue(std::unique_ptr<InsertionQueuer> ic) { Queue.push_back(std::move(ic)); };
+
     void create();
 
     void save_as(std::string path, std::string file_name) const;

src/Sketch/src/CircularArc.h

@@ -47,6 +47,8 @@ public:
 
     double length() const override { return radius() * arc_angle(); };
 
+    double length(double s0, double s1) const override { return std::abs(s1 - s0) * length(); };
+
     double radius() const { return Radius->value(); };
 
     void get_verticies(std::list<std::shared_ptr<Vertex const>> &v, MatchOrientation dir = MatchOrientation::Forward) const override {

src/Sketch/src/Curve.h

@@ -53,6 +53,8 @@ public: // interface
 
     virtual double length() const = 0; // length of segment between end points
 
+    virtual double length(double_t s0, double_t s1) const = 0; // length of segment between point s0 and s1
+
     virtual void get_verticies(std::list<std::shared_ptr<Vertex const>> &v, MatchOrientation dir = MatchOrientation::Forward) const = 0;
 
     virtual void replace_verticies(std::vector<std::shared_ptr<Vertex const>> const &oldv, std::vector<std::shared_ptr<Vertex const>> const &newv) = 0;

src/Sketch/src/LineSegment.h

@@ -25,6 +25,8 @@ public:
 
     double length() const override;
 
+    double length(double s0, double s1) const override { return std::abs(s1 - s0) * length(); };
+
     double area() const override { return 0.0; };
 
     double a(double s, bool orientation) const override;

test/Mesh/test_Mesh.cpp

@@ -346,13 +346,13 @@ TEST(Mesh, create_half_circle_domain) {
     m.save_as(SAVE_DIR, test_name + "_mesh");
 
     // Test number of verticies, edges, triangles
-    EXPECT_TRUE(m.size_points() == 6);
-    EXPECT_TRUE(m.size_edges() == 12);
-    EXPECT_TRUE(m.size_triangles() == 4);
+    EXPECT_EQ(m.size_points(), 10);
+    EXPECT_EQ(m.size_edges(), 24);
+    EXPECT_EQ(m.size_triangles(), 8);
 
-    EXPECT_TRUE(m.num_points() == 6);
-    EXPECT_TRUE(m.num_edges() == 9);
-    EXPECT_TRUE(m.num_triangles() == 4);
+    EXPECT_EQ(m.num_points(), 10);
+    EXPECT_EQ(m.num_edges(), 17);
+    EXPECT_EQ(m.num_triangles(), 8);
 
     // Test validity, optimality
     edges_are_valid(m);