Work on correcting fix for meshing errors still needs some cleanup

This may have been related to the delaunay optimality condition for triangles near the don't-care condition. There was a negative tolerance allowing sub-optimal elements to appear in the mesh which violated some assumptions and led to mesh divergence issues

src/Mesh/src/Mesh.cpp

@@ -374,7 +374,7 @@ bool Mesh::is_locally_optimal(size_t ei) const {
         double_t d2 = std::sqrt(v2x * v2x + v2y * v2y);
         double_t d3 = std::sqrt(v3x * v3x + v3y * v3y);
         double_t d4 = std::sqrt(v4x * v4x + v4y * v4y);
-        double_t tol = -std::sqrt(d1 * d2 * d3 * d4) * FLT_EPSILON;
+        //double_t tol = -std::sqrt(d1 * d2 * d3 * d4) * FLT_EPSILON * 0;
 
         double_t sina = v1x * v2y - v1y * v2x;
         double_t sinb = v3x * v4y - v3y * v4x;
@@ -382,7 +382,7 @@ bool Mesh::is_locally_optimal(size_t ei) const {
         double_t cosb = v3x * v4x + v3y * v4y;
         double_t cct = sina * cosb + cosa * sinb;
 
-        return cct >= tol; // TODO: Optional tol for testing purposes
+        return cct >= 0.0; // TODO: Optional tol for testing purposes
     }
 }
 
@@ -595,7 +595,7 @@ void Mesh::create() {
 
     insert_internal_boundaries();
 
-    get_triangles();
+    //get_triangles();
 }
 
 void Mesh::create_boundary_polygon() {
@@ -786,6 +786,7 @@ void Mesh::insert_internal_boundaries() {
         }
     }
 
+    /*
     // Attach edges to constraints by inserting interior curve midpoints until constraints are naturally satisfied
     std::list<size_t> queue;
     for (size_t i : interior_index) {
@@ -821,30 +822,6 @@ void Mesh::insert_internal_boundaries() {
             LocateTriangleResult result = locate_triangle(p0, ei);
 
             if (result != LocateTriangleResult::Point) {
-                /*
-                std::cerr << "Could not locate triangle referencing point, performing brute force search" << std::endl;
-                for (auto p : Points) {
-                    double_t d = hypot(p.X-p0.X,p.Y-p0.Y);
-                    //std::cout << d << std::endl;
-                }
-
-                bool point_found = false;
-                for (size_t i = 0; i != Edges.size(); ++i) {
-                    Point pe = point(Edges[i]);
-                    double_t d = hypot(pe.X-p0.X,pe.Y-p0.Y);
-                    if (d < DBL_MIN * 10.0) {
-                        std::cout << d << std::endl;
-                        ei = i;
-                        point_found = true;
-                        break;
-                    }
-                }
-
-                if(!point_found) {
-                    std::cerr << "Could not locate edge referencing point, aborting mesh" << std::endl;
-                    return;
-                }
-                 */
                 std::cerr << "Could not locate edge referencing point, aborting mesh" << std::endl;
                 return;
             }
@@ -885,12 +862,108 @@ void Mesh::insert_internal_boundaries() {
             }
         }
     }
+    */
+
+    // Insert all interior start and end points
+    double_t len_min = std::min(BoundaryConstraints[interior_index[0]]->curve()->length(),BoundaryConstraints[interior_index.back()]->curve()->length());
+    for (size_t i : interior_index) {
+        size_t N = (size_t)std::pow(2,std::ceil(std::log2(BoundaryConstraints[i]->curve()->length() / len_min / 2.0 + 1) + 1));
+        std::cout << N << std::endl;
+        for (size_t j = 0; j < N; ++j) {
+            double_t s = j * 1.0 / (N - 1);
+            Point p = BoundaryConstraints[i]->curve()->point(s);
+            LocateTriangleResult result = locate_triangle(p);
+            if (result == LocateTriangleResult::Interior) {
+                while (add_point_to_queue(p) == AddToQueueResult::Midpoint) {
+                    insert_from_queue();
+                }
+                insert_from_queue();
+            }
+        }
+        make_edges_optimal();
+        for (size_t j = 0; j < N - 1; ++j) {
+            double_t s0 = j * 1.0 / (N - 1);
+            Point p0 = BoundaryConstraints[i]->curve()->point(s0);
+
+            double_t s1 = (j+1) * 1.0 / (N - 1);
+            Point p1 = BoundaryConstraints[i]->curve()->point(s1);
+
+            size_t e{0};
+            locate_triangle(p1,e);
+            if (!find_attached(p0,e)) {
+                std::cout << "SOMETHING IS STILL WRONG" << std::endl;
+            }
+        }
+    }
+    return;
+
+    /*
+    while (queue.size() > 0) {
+        DartConstraint &dc = queue.back();
+        Point p0 = dc.curve()->point(dc.S0); // TODO: write Point Curve::point(double) and differentiate from Vertex Curve::vertex(double)
+        Point p1 = dc.curve()->point(dc.S1);
+
+        size_t ei = Edges.size() - 1;
+        LocateTriangleResult result = locate_triangle(p0, ei);
+
+        if (result == LocateTriangleResult::Exterior) {
+            std::cout << "Point is exterior to mesh, aborting" << std::endl;
+            std::cout << p0.X << "," << p0.Y << std::endl;
+            return;
+        } else if (result == LocateTriangleResult::Interior) {
+            std::cout << "Point is interior to mesh element, aborting" << std::endl;
+            std::cout << p0.X << "," << p0.Y << std::endl;
+            return;
+        }
+
+        if ((result == LocateTriangleResult::Point) && find_attached(p1, ei)) {
+            // if ei is an edge attaching p0 and p1, ei coincides with dc and is a constrained edge for BoundaryConstraint[i]
+            new_dart_constraint(dc.S0, dc.S1, dc.boundary_constraint());
+            DartConstraint &dcnew = DartConstraints.back();
+
+            Edge &e = Edges[ei];
+            e.Constraint = dcnew.Self;
+            e.Orientation = true;
+            dcnew.forward_dart(e.self());
+
+            Edge &et = Edges[e.Twin];
+            et.Constraint = dcnew.Self;
+            et.Orientation = false;
+            dcnew.reverse_dart(et.self());
+
+            queue.pop_back();
+
+            split_encroached_edges();
+        } else {
+            // if p0 and p1 are not connected by ei, insert the midpoint of dc
+            double_t s0 = dc.S0;
+            double_t s1 = dc.S1;
+            double_t sn = (s0 + s1) / 2.0;
+
+            Point const p = dc.curve()->point(sn);
+            AddToQueueResult atqr = add_point_to_queue(p);
+            if (atqr != AddToQueueResult::Duplicate) {
+                while (atqr == AddToQueueResult::Midpoint) {
+                    insert_from_queue();
+                    atqr = add_point_to_queue(p);
+                }
+                insert_from_queue();
+            } else {
+                std::cout << "Duplicate detected" << std::endl;
+            }
+
+            dc.S0 = sn;
+
+            queue.emplace_back(s0, sn, dc.boundary_constraint(), dc.self());
+        }
+    }
+     */
+    return;
 
     split_encroached_edges();
     for (auto bc : BoundaryConstraints) { // Length based refinement
         bc->refine(*this, 0.01);
     }
-
     split_encroached_edges();
 
     size_t np_old = Points.size();
@@ -909,6 +982,8 @@ void Mesh::insert_internal_boundaries() {
     }
 
     sort_constraints();
+
+    return;
 }
 
 void Mesh::make_edges_optimal() {
@@ -1083,8 +1158,8 @@ void Mesh::sort_constraints_by_length() {
      * Sorts BoundaryConstraints by length for internal boundary insertion
      */
 
-    //auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) { return (x->curve()->length()) < (y->curve()->length()); }; //ascending
-    auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) { return (x->curve()->length()) > (y->curve()->length()); }; //descending
+    auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) { return (x->curve()->length()) < (y->curve()->length()); }; // smallest to largest
+    //auto comp = [](std::shared_ptr<BoundaryConstraint> const &x, std::shared_ptr<BoundaryConstraint> const &y) { return (x->curve()->length()) > (y->curve()->length()); }; // largest to smallest
 
     std::sort(BoundaryConstraints.begin(), BoundaryConstraints.end(), comp);
 };
@@ -1312,6 +1387,80 @@ LocateTriangleResult Mesh::locate_triangle(Point const p, size_t &ei) const {
     double_t area12 = dx1p * dy2p - dx2p * dy1p;
     double_t area20 = dx2p * dy0p - dx0p * dy2p;
 
+    tol_a = FLT_EPSILON * (dx20 * dy01 - dy20 * dx01);
+    tol_l = FLT_EPSILON * sqrt(dx20 * dy01 - dy20 * dx01);
+
+    if (dist0 < tol_l) {
+        return LocateTriangleResult::Point;
+    } else if (dist1 < tol_l) {
+        ei = next(ei);
+        return LocateTriangleResult::Point;
+    } else if (dist2 < tol_l) {
+        ei = prev(ei);
+        return LocateTriangleResult::Point;
+    } else if (area01 >= 0.0 && area12 >= 0.0 && area20 >= 0.0) {
+        return LocateTriangleResult::Interior;
+    } else if (area01 < 0.0 && ei != twin(ei)) {
+        ei = twin(ei);
+        p2 = p1;
+        p1 = p0;
+        p0 = p2;
+
+        dx2p = dx1p;
+        dx1p = dx0p;
+        dx0p = dx2p;
+
+        dy2p = dy1p;
+        dy1p = dy0p;
+        dy0p = dy2p;
+
+        dx01 = -dx01;
+        dy01 = -dy01;
+
+        area01 = -area01;
+    } else if (area12 < 0.0 && next(ei) != twin(next(ei))) {
+        ei = twin(next(ei));
+        p0 = p2;
+
+        dx0p = dx2p;
+        dy0p = dy2p;
+
+        dx01 = -dx12;
+        dy01 = -dy12;
+
+        area01 = -area12;
+    } else if (area20 < 0.0 && prev(ei) != twin(prev(ei))) {
+        ei = twin(prev(ei));
+        p1 = p2;
+
+        dx1p = dx2p;
+        dy1p = dy2p;
+
+        dx01 = -dx20;
+        dy01 = -dy20;
+
+        area01 = -area20;
+    } else if (area01 < 0.0) {
+        return LocateTriangleResult::Exterior;
+    } else if (area12 < 0.0) {
+        ei = next(ei);
+        return LocateTriangleResult::Exterior;
+    } else if (area20 < 0.0) {
+        ei = prev(ei);
+        return LocateTriangleResult::Exterior;
+    } else {
+        std::cout << area01 << "," << area12 << "," << area20 << "," << tol_a << std::endl;
+        std::cout << xp << "," << yp << std::endl;
+        std::cout << p0.X << "," << p0.Y << std::endl;
+        std::cout << p1.X << "," << p1.Y << std::endl;
+        std::cout << p2.X << "," << p2.Y << std::endl;
+
+        std::cout << dist2 << "," << tol_l << std::endl;
+
+        throw std::runtime_error(std::string("Could not locate triangle containing point"));
+    }
+
+    /*
     if (dist0 < tol_l) {
         return LocateTriangleResult::Point;
     } else if (dist1 < tol_l) {
@@ -1374,8 +1523,10 @@ LocateTriangleResult Mesh::locate_triangle(Point const p, size_t &ei) const {
         ei = prev(ei);
         return LocateTriangleResult::Exterior;
     } else {
-        throw std::exception();
+        std::cerr << "Could not locate triangle containing point" << std::endl;
+        throw std::runtime_error(std::string("1"));
     }
+    */
 
     while (true) { // After first iteration, area01 > 0
         p2 = base(prev(ei));
@@ -1397,6 +1548,57 @@ LocateTriangleResult Mesh::locate_triangle(Point const p, size_t &ei) const {
         area12 = dx1p * dy2p - dx2p * dy1p;
         area20 = dx2p * dy0p - dx0p * dy2p;
 
+        if (dist2 < tol_l) {
+            ei = prev(ei);
+            return LocateTriangleResult::Point;
+        } else if (area12 >= 0.0 && area20 >= 0.0) {
+            //std::cout << area01 << "," << area12 << "," << area20 << std::endl;
+            return LocateTriangleResult::Interior;
+        } else if (area12 < 0.0 && next(ei) != twin(next(ei))) {
+                ei = twin(next(ei));
+                p0 = p2;
+
+                dx0p = dx2p;
+                dy0p = dy2p;
+
+                dx01 = -dx12;
+                dy01 = -dy12;
+
+                area01 = -area12;
+
+                continue;
+        } else if (area20 < 0.0 && prev(ei) != twin(prev(ei))) {
+                ei = twin(prev(ei));
+                p1 = p2;
+
+                dx1p = dx2p;
+                dy1p = dy2p;
+
+                dx01 = -dx20;
+                dy01 = -dy20;
+
+                area01 = -area20;
+
+                continue;
+        } else if (area12 < 0.0) {
+            ei = next(ei);
+            return LocateTriangleResult::Exterior;
+        } else if (area20 < 0.0) {
+            ei = prev(ei);
+            return LocateTriangleResult::Exterior;
+        } else {
+            std::cout << area01 << "," << area12 << "," << area20 << "," << tol_a << std::endl;
+            std::cout << xp << "," << yp << std::endl;
+            std::cout << p0.X << "," << p0.Y << std::endl;
+            std::cout << p1.X << "," << p1.Y << std::endl;
+            std::cout << p2.X << "," << p2.Y << std::endl;
+
+            std::cout << dist2 << "," << tol_l << std::endl;
+
+            throw std::runtime_error(std::string("Could not locate triangle containing point"));
+        }
+
+        /*
         if (dist2 < tol_l) {
             ei = prev(ei);
             return LocateTriangleResult::Point;
@@ -1435,8 +1637,17 @@ LocateTriangleResult Mesh::locate_triangle(Point const p, size_t &ei) const {
             ei = prev(ei);
             return LocateTriangleResult::Exterior;
         } else {
-            throw std::exception();
+            std::cout << area01 << "," << area12 << "," << area20 << "," << tol_a << std::endl;
+            std::cout << xp << "," << yp << std::endl;
+            std::cout << p0.X << "," << p0.Y << std::endl;
+            std::cout << p1.X << "," << p1.Y << std::endl;
+            std::cout << p2.X << "," << p2.Y << std::endl;
+
+            std::cout << dist2 << "," << tol_l << std::endl;
+
+            throw std::runtime_error(std::string("Could not locate triangle containing point"));
         }
+         */
     }
 }
 
@@ -1453,6 +1664,7 @@ InsertPointResult Mesh::insert_midpoint(size_t ei) {
     if (is_constrained(ei)) { // Constrained Edge
         c = DartConstraints.size();
 
+        size_t dci = Edges[ei].Constraint;
         DartConstraint &dc = DartConstraints[Edges[ei].Constraint];
         double_t s0 = dc.S0;
         double_t s1 = dc.S1;
@@ -1658,7 +1870,8 @@ AddToQueueResult Mesh::add_point_to_queue(Point vc, size_t ei) {
         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
+            std::cerr << "//TODO: No triangle found containing circumcenter, boundary edge was encroached by circumcenter, and the located edge was not a constrained edge" << std::endl;
+            throw std::runtime_error(std::string("3"));
         }
     }
 
@@ -1677,7 +1890,8 @@ AddToQueueResult Mesh::add_point_to_queue(Point vc, size_t ei) {
         return AddToQueueResult::Circumcenter;
 
     } else {
-        throw std::exception(); // TODO: No triangle could be found containing circumcenter and no boundary edge was encroached by circumcenter
+        std::cerr << "No triangle could be found containing circumcenter and no boundary edge was encroached by circumcenter" << std::endl;
+        throw std::runtime_error(std::string("4"));
     }
 }
 
@@ -1838,7 +2052,8 @@ std::shared_ptr<Contour const> Mesh::select_contour(Point const p) {
     }
 
     if ((*iter).edge() != e) { // Fail
-        throw std::exception();
+        std::cerr << "Failed to locate contour containing point" << std::endl;
+        throw std::runtime_error(std::string("4"));
     }
 
     return Contours[iter->contour()];

src/Sketch/src/Sketch.cpp

@@ -17,7 +17,7 @@
 
 #include <iostream>
 
-double Sketch::solve() {
+double_t Sketch::solve() {
     // TODO: Use sparse matrix representation for Jacobian, (this is lower priority since typical uses will not produce very large matrices)
     // TODO: Try extracting sub matrix and solving (Certain constraints are known to eliminate variables (Radius, Fixation), if the matrix can be permuted so that a submatrix on the block diagonal is lower-triangular, then a subsystem of equations can be solved independently of the remaining equations. Moreover, if the jacobian can be permuted into a block lower triangular form, subsets of equations can be solved blockwise)
     // TODO: Use variable_feature_size() and equation_feature_size()
@@ -28,19 +28,19 @@ double Sketch::solve() {
     Eigen::VectorXd lfs = Eigen::VectorXd::Zero(NumEquations);
 
     // Tolerances and iteration limits
-    double dbl_tol{std::sqrt(NumEquations + 1.0) * DBL_EPSILON};
-    double flt_tol{std::sqrt(NumEquations + 1.0) * FLT_EPSILON};
+    double_t dbl_tol{std::sqrt(NumEquations + 1.0) * DBL_EPSILON};
+    double_t flt_tol{std::sqrt(NumEquations + 1.0) * FLT_EPSILON};
 
-    size_t double_bits = (CHAR_BIT * sizeof(double));
+    size_t double_bits = (size_t)std::log2(1/DBL_EPSILON);
 
     // Initial residual calculation
     update_linearization(J, r);
 
     // Setup iteration
-    double res_norm{r.norm()};
-    double res_tol{characteristic_length() * dbl_tol};
+    double_t res_norm{r.norm()};
+    double_t res_tol{characteristic_length() * dbl_tol};
     long prev_rank(std::min(NumEquations, NumVariables));
-    for (size_t j = 0; res_norm > res_tol && j != double_bits; ++j) {
+    for (size_t j = 0; res_norm > res_tol && j < double_bits; ++j) {
         Eigen::FullPivHouseholderQR<Eigen::MatrixXd> QR = J.fullPivHouseholderQr();
         Eigen::VectorXd delta = QR.solve(r);
 
@@ -57,8 +57,8 @@ double Sketch::solve() {
         update_variables(delta);
         update_linearization(J, r);
 
-        double iter_norm{r.norm()};
-        double iter_tol{res_tol};
+        double_t iter_norm{r.norm()};
+        double_t iter_tol{res_tol};
         if (iter_norm > res_norm) { // Line search
             delta *= -0.5;
 
@@ -74,7 +74,8 @@ double Sketch::solve() {
                     line_search_iter += 1;
 
                     if (line_search_iter == double_bits) {
-                        std::cerr << "Residual norm was not reduced for any step-size greater than 2^-(CHAR_BIT * sizeof(double))" << std::endl;
+                        std::cerr << "Residual norm was not reduced for any step-size greater than DBL_EPSILON" << std::endl;
+                        break;
                     }
                 }
             }
@@ -108,7 +109,7 @@ void Sketch::update_variables(Eigen::VectorXd &delta) const {
     }
 }
 
-void Sketch::perturb(Eigen::VectorXd &delta, double scale) const {
+void Sketch::perturb(Eigen::VectorXd &delta, double_t scale) const {
     // Randomly perturb the delta vector to prevent stagnation at inflection point
     std::random_device rdev;
     std::mt19937 rgen(rdev());
@@ -119,8 +120,8 @@ void Sketch::perturb(Eigen::VectorXd &delta, double scale) const {
     }
 }
 
-double Sketch::characteristic_length() const {
-    double scale{0.0};
+double_t Sketch::characteristic_length() const {
+    double_t scale{0.0};
 
     if (Curves.size() > 0) {
         for (size_t i = 0; i != Curves.size(); ++i) {

test/Model Templates/test_SynchronousReluctanceMachine.cpp

@@ -44,7 +44,7 @@ TEST_F(Synchronous_Reluctance_Machine, geometry_test) {
     srr.add_to_sketch(sketch,origin);
     dws.add_to_sketch(sketch,origin);
 
-    CylindricalAirgap<1> cag{sketch,origin,srr.OuterRadius,0.0,dws.InnerRadius,0.0,360.0/8.0};
+    CylindricalAirgap<2> cag{sketch,origin,srr.OuterRadius,0.0,dws.InnerRadius,0.0,360.0/8.0};
 
     EXPECT_LE(sketch.solve(), 100e-3 * FLT_EPSILON);
     EXPECT_TRUE(sketch.build());
@@ -58,9 +58,9 @@ TEST_F(Synchronous_Reluctance_Machine, geometry_test) {
     EXPECT_TRUE(mesh.edges_are_valid());
     EXPECT_TRUE(mesh.edges_are_optimal());
 
-    mesh.MinimumElementQuality = 0.5;
+    //mesh.MinimumElementQuality = 0.5;
 
-    EXPECT_TRUE(mesh.refine());
+    //EXPECT_TRUE(mesh.refine());
 
     mesh.save_as(SAVE_DIR,"Synchronous_Reluctance_Machine");
 }
\ No newline at end of file