FIx problem with Sketch Contour/Boundary parsing related to supremum calculations

.idea/Oersted.iml

@@ -74,6 +74,7 @@
       <sourceFolder url="file://$MODULE_DIR$/test/Mesh/test_Mesh.cpp" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/test/Mesh/test_Mesh.hpp" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/test/CMakeLists.txt" isTestSource="false" />
+      <sourceFolder url="file://$MODULE_DIR$/test/Sketch/test_Curve.cpp" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/test/Sketch/test_LineSegment.cpp" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/test/Sketch/test_Pattern.cpp" isTestSource="false" />
       <sourceFolder url="file://$MODULE_DIR$/test/Sketch/test_Sketch.cpp" isTestSource="false" />

.idea/dictionaries/jpries.xml

@@ -2,6 +2,7 @@
   <dictionary name="jpries">
     <words>
       <w>eigen</w>
+      <w>supremum</w>
       <w>tangency</w>
       <w>verticies</w>
     </words>

src/Sketch/src/CircularArc.cpp

@@ -97,7 +97,7 @@ double CircularArc::arc_angle() const {
 Vertex CircularArc::point(double s) const {
     double a = s_to_a(s);
 
-    return Vertex{Center->x() + radius() * cos(a), Center->y() + radius() * sin(a)};
+    return Vertex{center()->x() + radius() * cos(a), center()->y() + radius() * sin(a)};
 }
 
 Vertex CircularArc::tangent(double s, bool orientation) const {
@@ -149,15 +149,41 @@ double CircularArc::da(double s, bool orientation) const {
     }
 }
 
-double CircularArc::supremum() const {
-    double sup = std::fmax(Start->hypot(), End->hypot());
-    double s = a_to_s(Center->atan());
+std::pair<double, double> CircularArc::supremum() const {
+    double x = start()->x();
+    double y = start()->y();
+    double sup = sqrt(x * x + y * y);
+    double par = 0.0;
+
+    double xx = end()->x();
+    double yy = end()->y();
+    double val = sqrt(xx * xx + yy * yy);
+    if (val > sup) {
+        x = xx;
+        y = yy;
+        sup = val;
+        par = 1.0;
+    }
 
-    if (s > 0 && s < 1) {
-        sup = std::fmax(sup, point(s).hypot());
+    double s = a_to_s(center()->atan());
+    if (s > 0.0 && s < 1.0) {
+        Vertex v = point(s);
+        xx = v.x();
+        yy = v.y();
+        val = sqrt(xx * xx + yy * yy);
+
+        if(val > sup) {
+            x = xx;
+            y = yy;
+            sup = val;
+            par = s;
+        }
     }
 
-    return std::move(sup);
+    double ang = s_to_a(par);
+    double cross = abs(x * cos(ang) + y * sin(ang)) / sup; // cross product of vector from origin to point and tangent vector
+
+    return std::pair<double,double>(sup, cross);
 }
 
 bool CircularArc::on_manifold(const double x, const double y) const {

src/Sketch/src/CircularArc.h

@@ -58,7 +58,7 @@ public:
 
     double da(double s, bool orientation) const override;
 
-    double supremum() const override;
+    std::pair<double, double> supremum() const override;
 
     // Curve-Vertex Comparison
     using Curve::on_manifold;

src/Sketch/src/Constellation.cpp

@@ -26,14 +26,14 @@ bool Constellation::twin(std::list<Star>::iterator &s_out, std::list<Branch>::it
     return false;
 }
 
-void Constellation::supremum(std::list<Star>::iterator &s_out, std::list<Branch>::iterator &b_out) {
-    double sup = 0.0;
-    double ang;
+void Constellation::supremum(list<Star>::iterator &s_out, list<Branch>::iterator &b_out) {
+    pair<double,double> sup{0.0,0.0};
+    double ang{0.0};
 
     for (auto s = Stars.begin(); s != Stars.end(); ++s) {
         for (auto b = s->begin(); b != s->end(); ++b) {
-            double sup_ij = b->Path->supremum();
-            if (sup < sup_ij || (sup == sup_ij && b->Angle < ang)) {
+            std::pair<double,double> sup_ij = b->Path->supremum();
+            if ((sup_ij > sup) || (sup == sup_ij && b->Angle < ang)) {
                 sup = sup_ij;
                 ang = b->Angle;
 
@@ -76,14 +76,14 @@ void Constellation::pop(const Curve *c) {
 }
 
 bool Constellation::boundary(Contour *c) {
-    std::vector<const Curve *> curves;
-    std::vector<bool> orientation;
+    vector<const Curve *> curves;
+    vector<bool> orientation;
 
     // Base of induction
     auto s{Stars.begin()};
     auto b{s->begin()};
 
-    supremum(s, b);
+    supremum(s,b);
 
     double angle{0.0};
     b = s->prev(b);
@@ -122,8 +122,8 @@ bool Constellation::boundary(Contour *c) {
 }
 
 bool Constellation::contours(std::vector<Contour *> &contours) {
-    std::vector<const Curve *> contour_curves;
-    std::vector<bool> orientation;
+    vector<const Curve *> contour_curves;
+    vector<bool> orientation;
 
     while (size() > 0) {
         bool success = find_closed_contour(contour_curves, orientation);
@@ -145,7 +145,7 @@ bool Constellation::find_closed_contour(std::vector<const Curve *> &curves, std:
     auto s{Stars.begin()};
     auto b{s->begin()};
 
-    supremum(s, b);
+    supremum(s,b);
 
     double angle{0.0};
     b = s->next(b);

src/Sketch/src/Constellation.h

@@ -12,20 +12,20 @@ public:
 
     size_t size() { return Stars.size(); };
 
-    bool contours(std::vector<Contour *> &contours);
+    bool contours(vector<Contour *> &contours);
 
     bool boundary(Contour *c);
 
 private:
-    std::list<Star> Stars;
+    list<Star> Stars;
 
     void pop(const Curve *c = nullptr);
 
-    bool twin(std::list<Star>::iterator &s_out, std::list<Branch>::iterator &b_out);
+    bool twin(list<Star>::iterator &s_out, list<Branch>::iterator &b_out);
 
-    void supremum(std::list<Star>::iterator &s_out, std::list<Branch>::iterator &b_out);
+    void supremum(list<Star>::iterator &s_out, list<Branch>::iterator &b_out);
 
-    bool find_closed_contour(std::vector<const Curve *> &curves, std::vector<bool> &orientation);
+    bool find_closed_contour(vector<const Curve *> &curves, vector<bool> &orientation);
 };
 
 #endif //OERSTED_CONSTELLATION_H
\ No newline at end of file

src/Sketch/src/Curve.h

@@ -33,10 +33,10 @@ public:
     virtual double length() const = 0; // length of segment between end points
     virtual double area() const = 0;    // area of segment between end points
 
-    virtual double a(double s, bool orientation) const = 0;        // tangent angle
+    virtual double a(double s, bool orientation) const = 0;     // tangent angle
     virtual double da(double s, bool orientation) const = 0;    // derivative of tangent angle per unit length
 
-    virtual double supremum() const = 0;    // maximum length of vector between origin and point on curve
+    virtual std::pair<double, double> supremum() const = 0;     // maximum length of vector between origin and point on curve
 
     // Curve-Vertex Comparison
     virtual bool on_manifold(const Vertex *v) const final; // true if vertex is on manifold defined by curve

src/Sketch/src/LineSegment.cpp

@@ -32,8 +32,31 @@ double LineSegment::a(double s, bool orientation) const {
     }
 };
 
-double LineSegment::supremum() const {
-    return std::fmax(start()->hypot(), end()->hypot());
+std::pair<double, double> LineSegment::supremum() const {
+    double xs = start()->x();
+    double ys = start()->y();
+    double ls = sqrt(xs * xs + ys * ys);
+
+    double xe = end()->x();
+    double ye = end()->y();
+    double le = sqrt(xe * xe + ye * ye);
+
+    double dx = xe - xs;
+    double dy = ye - ys;
+    double dl = sqrt(dx * dx + dy * dy);
+    dx /= dl;
+    dy /= dl;
+
+    double sup, cross;
+    if (ls >= le) {
+        sup = ls;
+        cross = abs(xs * dy - ys * dx) / ls;
+    } else {
+        sup = le;
+        cross = abs(xe * dy - ye * dx) / le;
+    }
+
+    return std::pair<double,double>(sup,cross);
 };
 
 bool LineSegment::on_manifold(const double x, const double y) const {

src/Sketch/src/LineSegment.h

@@ -54,7 +54,7 @@ public:
 
     double da(double s, bool orientation) const override { return 0.0; };
 
-    double supremum() const override;
+    std::pair<double, double> supremum() const override;
 
     // Curve-Vertex Comparison
     using Curve::on_manifold;

src/Sketch/src/Sketch.cpp

@@ -50,6 +50,7 @@ void Sketch::solve() {
 bool Sketch::build() {
     Constellation c = Constellation(this);
 
+
     bool success = c.boundary(Boundary);
 
     Contours.resize(0);

src/Sketch/src/Sketch.h

@@ -14,6 +14,8 @@
 
 #include <boost/filesystem.hpp>
 
+using namespace std;
+
 class Sketch;
 
 // Sketch Parameter

test/CMakeLists.txt

@@ -8,9 +8,12 @@ set(SOURCE_FILES
         Sketch/test_Vertex.cpp
         Sketch/test_LineSegment.cpp
         Sketch/test_CircularArc.cpp
+        Sketch/test_Curve.cpp
+
         Sketch/test_Star.cpp
         Sketch/test_Contour.cpp
         Sketch/test_Sketch.cpp
+
         Sketch/test_Constraint.cpp
         Sketch/test_Pattern.cpp
 

test/Mesh/test_Mesh.cpp

@@ -704,7 +704,8 @@ TEST(MESH__LOCATE_TRIANGLE, TRIANGULAR_DOMAIN) {
     LineSegment &l2 = s.new_element<LineSegment>(v2, v0);
 
     s.solve();
-    s.build();
+    EXPECT_TRUE(s.build());
+
     /*
     std::vector<const Curve*> cc{ &l0,&l1,&l2 };
     Contour cont{ cc };

test/Sketch/test_Constraint.cpp

@@ -432,7 +432,7 @@ TEST(Constraint, Distance_LineSegment_midpoint_intersection) {
 
     double dot = abs(dx0 * dx1 + dy0 * dy1);
 
-    EXPECT_NEAR(1.0, dot, TOL); // #TODO: Fails because lines intersect at midpoint
+    EXPECT_NEAR(1.0, dot, TOL); // TODO: Fails because lines intersect at midpoint
 }
 
 TEST(Constraint, Distance_CircularArc_exterior) {

test/Sketch/test_Contour.cpp

@@ -136,8 +136,8 @@ TEST(CONTOUR, IMPLICIT_SELF_INTERSECTION_FAILURE) {
     EXPECT_ANY_THROW(Contour cont{c}); // Construction should fail since contour is self intersecting
 
     /*
-        #TODO: Use nurbs representation to approximate potential intersection point.
-        #TODO: Use newton's method to solve intersection problem with the generated initial guesses.
+        TODO: Use nurbs representation to approximate potential intersection point.
+        TODO: Use newton's method to solve intersection problem with the generated initial guesses.
     */
 }
 

test/Sketch/test_Curve.cpp

+#include "test_Sketch.hpp"
+
+TEST(Curve, Supremum) {
+    Sketch s;
+
+    Vertex &origin = s.new_element<Vertex>(0.0, 0.0);
+    Vertex &v0 = s.new_element<Vertex>(1.0, 0.0);
+    Vertex &v1 = s.new_element<Vertex>(4.0, 0.0);
+    Vertex &v2 = s.new_element<Vertex>(M_SQRT2, M_SQRT2);
+    Vertex &v3 = s.new_element<Vertex>(M_SQRT1_2, M_SQRT1_2);
+
+    LineSegment &l0 = s.new_element<LineSegment>(v0, v1);
+    LineSegment &l1 = s.new_element<LineSegment>(v3, v2);
+    CircularArc &c0 = s.new_element<CircularArc>(v1, v2, origin, 4.0);
+    CircularArc &c1 = s.new_element<CircularArc>(v0, v3, origin, 1.0);
+
+    Radius &r0 = s.new_element<Radius>(c0, 4.0);
+    Radius &r1 = s.new_element<Radius>(c1, 1.0);
+    Horizontal &h0 = s.new_element<Horizontal>(l0);
+    Angle &a0 = s.new_element<Angle>(l0, l1, 45.0);
+    Coincident<LineSegment> &coin0 = s.new_element<Coincident<LineSegment>>(origin, l0);
+    Coincident<LineSegment> &coin1 = s.new_element<Coincident<LineSegment>>(origin, l1);
+
+    std::pair<double,double> sc0 = c0.supremum();
+    std::pair<double,double> sl0 = l0.supremum();
+
+    EXPECT_GT(sc0, sl0);
+
+    s.solve();
+    s.build();
+}
\ No newline at end of file

test/Sketch/test_Star.cpp

@@ -210,8 +210,8 @@ TEST(STAR, FIND_CLOSED_CONTOUR_0) {
     // Sketch internal contour parsing
     {
         sketch.solve();
-        sketch.build();
-        EXPECT_TRUE(sketch.size_contours() == 1);
+        EXPECT_TRUE(sketch.build());
+        EXPECT_EQ(sketch.size_contours(), 1);
         EXPECT_TRUE(*sketch.contour(0) == *sketch.boundary());
 
         const Contour *contour = sketch.contour(0);
@@ -296,22 +296,22 @@ TEST(STAR, FIND_CLOSED_CONTOUR_2) {
     // Sketch internal contour parsing
     {
         sketch.solve();
-        sketch.build();
+        EXPECT_TRUE(sketch.build());
 
-        EXPECT_TRUE(sketch.size_contours() == 2);
+        EXPECT_EQ(sketch.size_contours(), 2);
         EXPECT_FALSE(sketch.contour(0) == sketch.boundary());
         EXPECT_FALSE(sketch.contour(1) == sketch.boundary());
 
         const Contour *contour = sketch.contour(1);
 
-        EXPECT_TRUE(contour->size() == 3);
+        EXPECT_EQ(contour->size(), 3);
         EXPECT_TRUE(&l0 == contour->curve(0) || &l0 == contour->curve(1) || &l0 == contour->curve(2));
         EXPECT_TRUE(&l1 == contour->curve(0) || &l1 == contour->curve(1) || &l1 == contour->curve(2));
         EXPECT_TRUE(&c0 == contour->curve(0) || &c0 == contour->curve(1) || &c0 == contour->curve(2));
 
         contour = sketch.contour(0);
 
-        EXPECT_TRUE(contour->size() == 3);
+        EXPECT_EQ(contour->size(), 3);
         EXPECT_TRUE(&l3 == contour->curve(0) || &l3 == contour->curve(1) || &l3 == contour->curve(2));
         EXPECT_TRUE(&l4 == contour->curve(0) || &l4 == contour->curve(1) || &l4 == contour->curve(2));
         EXPECT_TRUE(&c0 == contour->curve(0) || &c0 == contour->curve(1) || &c0 == contour->curve(2));
@@ -331,15 +331,11 @@ TEST(STAR, FIND_CLOSED_CONTOUR_2) {
     {
         const Contour *boundary = sketch.boundary();
 
-        EXPECT_TRUE(boundary->size() == 4);
-        EXPECT_TRUE(&l0 == boundary->curve(0) || &l0 == boundary->curve(1) || &l0 == boundary->curve(2) ||
-                    &l0 == boundary->curve(3));
-        EXPECT_TRUE(&l1 == boundary->curve(0) || &l1 == boundary->curve(1) || &l1 == boundary->curve(2) ||
-                    &l1 == boundary->curve(3));
-        EXPECT_TRUE(&l3 == boundary->curve(0) || &l3 == boundary->curve(1) || &l3 == boundary->curve(2) ||
-                    &l3 == boundary->curve(3));
-        EXPECT_TRUE(&l4 == boundary->curve(0) || &l4 == boundary->curve(1) || &l4 == boundary->curve(2) ||
-                    &l4 == boundary->curve(3));
+        EXPECT_EQ(boundary->size(), 4);
+        EXPECT_TRUE(&l0 == boundary->curve(0) || &l0 == boundary->curve(1) || &l0 == boundary->curve(2) || &l0 == boundary->curve(3));
+        EXPECT_TRUE(&l1 == boundary->curve(0) || &l1 == boundary->curve(1) || &l1 == boundary->curve(2) || &l1 == boundary->curve(3));
+        EXPECT_TRUE(&l3 == boundary->curve(0) || &l3 == boundary->curve(1) || &l3 == boundary->curve(2) || &l3 == boundary->curve(3));
+        EXPECT_TRUE(&l4 == boundary->curve(0) || &l4 == boundary->curve(1) || &l4 == boundary->curve(2) || &l4 == boundary->curve(3));
     }
 }
 

test/UseCases/test_Stator.cpp

@@ -35,13 +35,13 @@ TEST(STATOR, 0) {
 
     //std::vector<const Curve*> mv{ &l0,&l1,&l2,&l3,&c0,&c1,&c2,&c3 };
     std::vector<const Curve*> mv{&l0, &l3, &c0, &c1, &c2, &c3};
-    // l1.ForConstruction = true; //TODO: Fixes segfault (#1) failure but should not be necessary, problem is related to MirrorCopy mirror line l1
     MirrorCopy &m0 = sketch.new_element<MirrorCopy>(mv, &l1);
 
     sketch.solve();
 
     sketch.save_as<SaveMethod::Rasterize>(SAVE_DIR, "stator0d0.csv");
 
+    EXPECT_TRUE(sketch.build());
     EXPECT_EQ(sketch.boundary()->size(), 8);
 
     Fixation &f = sketch.new_element<Fixation>(origin);
@@ -66,17 +66,15 @@ TEST(STATOR, 0) {
 
     sketch.save_as<SaveMethod::Rasterize>(SAVE_DIR, "stator0d1.csv");
 
+    EXPECT_TRUE(sketch.build());
     EXPECT_EQ(sketch.boundary()->size(), 8);
 
-    sketch.build();
-
     Mesh mesh{sketch};
 
     mesh.MinimumElementQuality = M_SQRT1_2;
     mesh.MaximumElementSize = 2.5;
     mesh.MinimumElementSize = 0.25;
 
-    /* TODO: Segfault (#1)
     mesh.create();
 
     mesh.save_as(SAVE_DIR, "stator_0.csv");
@@ -84,5 +82,4 @@ TEST(STATOR, 0) {
     mesh.refine();
 
     mesh.save_as(SAVE_DIR, "stator_0_refine.csv");
-    */
 }
\ No newline at end of file