Adds Sketch methods for selecting periodic and radial boundaries

src/Sketch/CMakeLists.txt

@@ -16,6 +16,8 @@ set(SOURCE_FILES
             ./src/LineSegment.h ./src/LineSegment.cpp
             ./src/CircularArc.h ./src/CircularArc.cpp
 
+        ./src/ContinuousBoundaryPair.h ./src/ContinuousBoundaryPair.cpp
+
         ./src/Constraint.h  ./src/Constraint.cpp
             ./src/Angle.h           ./src/Angle.cpp
             ./src/Coincident.h    ./src/Coincident.cpp

src/Sketch/src/ContinuousBoundaryPair.cpp

+#include "ContinuousBoundaryPair.h"

src/Sketch/src/ContinuousBoundaryPair.h

+#ifndef OERSTED_CONTINUOUSBOUNDARYPAIR_H
+#define OERSTED_CONTINUOUSBOUNDARYPAIR_H
+
+#include "Curve.h"
+
+class ContinuousBoundaryPair {
+public:
+    ContinuousBoundaryPair(std::shared_ptr<Curve const> c0, std::shared_ptr<Curve const> c1, bool dir) : Curve0{c0}, Curve1{c1}, Orientation{dir} {};
+
+    std::shared_ptr<Curve const> Curve0;
+    std::shared_ptr<Curve const> Curve1;
+
+    bool Orientation;
+};
+
+#endif //OERSTED_CONTINUOUSBOUNDARYPAIR_H

src/Sketch/src/Contour.h

@@ -36,6 +36,8 @@ public:
 
     std::shared_ptr<Curve const> curve(size_t i) const { return Curves[i]; };
 
+    std::vector<std::shared_ptr<Curve const>> curves() const { return Curves; };
+
 private:
     std::vector<std::shared_ptr<Curve const>> Curves;
 

src/Sketch/src/Sketch.cpp

 #include <random>
+#include <utility>
 
 #include "Sketch.h"
 #include "Vertex.h"
+#include "CircularArc.h"
+#include "LineSegment.h"
 #include "Curve.h"
 #include "Constraint.h"
 #include "Pattern.h"
@@ -198,4 +201,63 @@ void Sketch::save_as<SaveMethod::Rasterize>(std::string path, std::string file_n
     }
 
     fs.close();
+}
+
+std::vector<ContinuousBoundaryPair> Sketch::select_periodic_boundary_pairs(std::shared_ptr<Vertex> v0, double_t angle) const {
+    std::vector<ContinuousBoundaryPair> cbp;
+
+    std::vector<std::shared_ptr<Curve const>> bc = Boundary->curves();
+
+    while (bc.begin() != bc.end()) {
+        auto c0 = bc.begin();
+        auto c1 = --bc.end();
+
+        MatchOrientation m{MatchOrientation::None};
+        while (c0 != c1) {
+            auto cc0 = (*c0).get();
+            auto cc1 = (*c1).get();
+            m = (*c1)->is_identical(*c0, v0, angle);
+
+            if (m == MatchOrientation::None) {
+                m = (*c0)->is_identical(*c1, v0, angle);
+                if (m != MatchOrientation::None) {
+                    std::swap(c0, c1);
+                }
+            }
+
+            if (m == MatchOrientation::Forward) {
+                cbp.emplace_back(*c0, *c1, true);
+                break;
+            } else if (m == MatchOrientation::Reverse) {
+                cbp.emplace_back(*c0, *c1, false);
+                break;
+            } else {
+                ++c0;
+            }
+        }
+
+        if (m != MatchOrientation::None) { // if match found, erase c1
+            bc.erase(c0);
+            bc.erase(c1);
+        } else {
+            bc.erase(c1);
+        }
+    }
+
+    return cbp;
+}
+
+std::vector<std::shared_ptr<Curve const>> Sketch::select_radial_boundary(std::shared_ptr<Vertex> v0, double_t radius) const {
+    std::vector<std::shared_ptr<Curve const>> rb;
+
+    std::vector<std::shared_ptr<Curve const>> bc = Boundary->curves();
+
+    auto test_circle = std::make_shared<CircularArc>(std::make_shared<Vertex>(), std::make_shared<Vertex>(), v0, radius);
+    for (std::shared_ptr<Curve const> c : bc) {
+        if (test_circle->is_coincident(c)) {
+            rb.push_back(c);
+        }
+    }
+
+    return rb;
 }
\ No newline at end of file

src/Sketch/src/Sketch.h

@@ -14,6 +14,8 @@
 
 #include "Eigen"
 
+#include "ContinuousBoundaryPair.h"
+
 class Variable;
 class Vertex;
 class Curve;
@@ -81,6 +83,10 @@ public:
 
     std::vector<std::shared_ptr<Curve const>> const curves() { return std::vector<std::shared_ptr<Curve const>>{Curves.begin(),Curves.end()};};
 
+    std::vector<ContinuousBoundaryPair> select_periodic_boundary_pairs(std::shared_ptr<Vertex> v0, double_t angle) const;
+
+    std::vector<std::shared_ptr<Curve const>> select_radial_boundary(std::shared_ptr<Vertex> v0, double_t radius) const;
+
 private:
     std::vector<std::shared_ptr<Variable>> Variables;
 

test/LibraryIntegration/test_Mesh_to_FEM.cpp

 #include "test_Library_Integration.hpp"
 
-TEST(Library_Integration_Circle, Uniform_Current_Density) {
-    std::string save_dir = "./test/output/Integration/";
+std::string SAVE_DIR = "./test/output/Integration/";
 
+TEST(Library_Integration, Full_Circle_Uniform_Current_Density) {
     // Create Sketch
     Sketch sk;
 
-    auto v0 = sk.new_element<Vertex>(0.0,0.0);
-    auto v1 = sk.new_element<Vertex>(1.0,0.0);
+    auto v0 = sk.new_element<Vertex>(0.0, 0.0);
+    auto v1 = sk.new_element<Vertex>(1.0, 0.0);
     auto v2 = sk.new_element<Vertex>(1.0 * std::cos(M_PI * 2.0 / 3.0), 1.0 * std::sin(M_PI * 2.0 / 3.0));
     auto v3 = sk.new_element<Vertex>(1.0 * std::cos(-M_PI * 2.0 / 3.0), 1.0 * std::sin(-M_PI * 2.0 / 3.0));
 
-    auto c0 = sk.new_element<CircularArc>(v1,v2,v0,1.0);
-    auto c1 = sk.new_element<CircularArc>(v2,v3,v0,1.0);
-    auto c2 = sk.new_element<CircularArc>(v3,v1,v0,1.0);
+    auto c0 = sk.new_element<CircularArc>(v1, v2, v0, 1.0);
+    auto c1 = sk.new_element<CircularArc>(v2, v3, v0, 1.0);
+    auto c2 = sk.new_element<CircularArc>(v3, v1, v0, 1.0);
 
     auto f0 = sk.new_element<Fixation>(v0);
     auto f1 = sk.new_element<Fixation>(v1);
@@ -26,9 +26,9 @@ TEST(Library_Integration_Circle, Uniform_Current_Density) {
     bool result = sk.build();
     EXPECT_TRUE(result);
 
-    EXPECT_EQ(sk.size_contours(),1);
+    EXPECT_EQ(sk.size_contours(), 1);
 
-    sk.save_as<SaveMethod::Rasterize>(save_dir, std::string("circle_sketch"));
+    sk.save_as<SaveMethod::Rasterize>(SAVE_DIR, std::string("circle_sketch"));
 
     // Create Refineable Mesh
     Mesh me{sk};
@@ -40,21 +40,21 @@ TEST(Library_Integration_Circle, Uniform_Current_Density) {
 
     me.refine();
 
-    me.save_as(save_dir, std::string("circle_mesh"));
+    me.save_as(SAVE_DIR, std::string("circle_mesh"));
 
     // Convert to FiniteElementMesh
     FiniteElementMesh<2, 1> fem{me};
 
-    for (std::shared_ptr<Boundary<2>> const &b : fem.boundaries() ){
-        double_t a0{-2*M_PI};
-        for(size_t i : b->nodes()) {
+    for (std::shared_ptr<Boundary<2>> const &b : fem.boundaries()) {
+        double_t a0{-2 * M_PI};
+        for (size_t i : b->nodes()) {
             XY const &p = fem.node(i);
 
             // Test radius of boundary
             EXPECT_NEAR(1.0, std::hypot(p.x(), p.y()), FLT_EPSILON);
 
             // Test ordering of boundary nodes
-            double_t a1{std::atan2(p.y(),p.x())};
+            double_t a1{std::atan2(p.y(), p.x())};
             if (a1 < 0) {
                 a1 += 2 * M_PI;
             }
@@ -66,7 +66,7 @@ TEST(Library_Integration_Circle, Uniform_Current_Density) {
     // Create magnetostatic physics
     Magnetostatic<2, 1, 1, FieldVariable::A> msph{fem};
 
-    msph.add_current_density([](double t) { return 1.0 / (2.0*M_PI*1e-7); }, {fem.region(0)});
+    msph.add_current_density([](double t) { return 1.0 / (2.0 * M_PI * 1e-7); }, {fem.region(0)});
 
     msph.add_magnetic_insulation({fem.boundary(0), fem.boundary(1), fem.boundary(2)});
 
@@ -126,13 +126,95 @@ TEST(Library_Integration_Circle, Uniform_Current_Density) {
     for (size_t i = 0; i != qp.size(); ++i) {
         double_t r = std::hypot(qp[i][0].x(), qp[i][0].y());
 
-        EXPECT_NEAR(Bmag(i), r, 0.1);
+        EXPECT_NEAR(Bmag(i), r, 1.0 * 0.05);
 
         double_t a = std::atan2(qp[i][0].y(), qp[i][0].x()) * 180 / M_PI + 90.0;
 
         if (a > 180.0) { a -= 360.0; }
-        if (a - Bang(i) > 180.0)  { a -= 360.0; }
+        if (a - Bang(i) > 180.0) { a -= 360.0; }
         if (a - Bang(i) < -180.0) { a += 360.0; }
-        EXPECT_NEAR(Bang(i), a, 360 * 0.1);
+        EXPECT_NEAR(Bang(i), a, 360 * 0.05);
     }
+}
+
+TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) {
+    // Create Sketch
+    Sketch sk;
+
+    auto v0 = sk.new_element<Vertex>(0.0, 0.0);
+    auto v1 = sk.new_element<Vertex>(1.0, 0.0);
+    auto v2 = sk.new_element<Vertex>(2.0, 0.0);
+    auto v3 = sk.new_element<Vertex>(M_SQRT1_2, M_SQRT1_2);
+    auto v4 = sk.new_element<Vertex>(M_SQRT2, M_SQRT2);
+
+    auto f0 = sk.new_element<Fixation>(v0);
+
+    auto l01 = sk.new_element<LineSegment>(v0, v1);
+    auto l12 = sk.new_element<LineSegment>(v1, v2);
+
+    auto l04 = sk.new_element<LineSegment>(v0, v4, true);
+
+    auto h01 = sk.new_element<Horizontal>(l01);
+    auto h12 = sk.new_element<Horizontal>(l12);
+
+    auto a0103 = sk.new_element<Angle>(l01, l04, 45.0);
+
+    auto c013 = sk.new_element<CircularArc>(v1, v3, v0, 1.0);
+    auto c024 = sk.new_element<CircularArc>(v2, v4, v0, 2.0);
+
+    auto r013 = sk.new_element<Radius>(c013, 1.0);
+    auto r024 = sk.new_element<Radius>(c024, 2.0);
+
+    std::vector<std::shared_ptr<Curve const>> vec{l01, l12, c013, c024};
+    auto m04 = sk.new_element<MirrorCopy>(vec, l04);
+
+    double residual = sk.solve();
+    EXPECT_LE(residual, FLT_EPSILON);
+
+    bool result = sk.build();
+    ASSERT_TRUE(result);
+
+    sk.save_as<SaveMethod::Rasterize>(SAVE_DIR, "quarter_circle_mirror_copy_sketch");
+
+    auto periodic_boundary = sk.select_periodic_boundary_pairs(v0, 90.0);
+    {
+        for (auto &bp : periodic_boundary) {
+            if (bp.Curve0.get() == l01.get()) {
+                EXPECT_EQ(bp.Curve1->is_identical(l01, v0, 90.0), MatchOrientation::Reverse);
+            } else if (bp.Curve0.get() == l12.get()) {
+                EXPECT_EQ(bp.Curve1->is_identical(l12, v0, 90.0), MatchOrientation::Reverse);
+            } else {
+                GTEST_NONFATAL_FAILURE_("No matching boundary found");
+            }
+        }
+    }
+
+    auto radial_boundary = sk.select_radial_boundary(v0, 2.0);
+    {
+        for (auto &c : radial_boundary) {
+            auto cc = std::dynamic_pointer_cast<CircularArc const>(c);
+            if (cc) {
+                EXPECT_NEAR(cc->radius(), 2.0, FLT_EPSILON);
+                EXPECT_EQ(cc->center().get(), v0.get());
+            } else {
+                GTEST_NONFATAL_FAILURE_("dynamic_cast of Curve to CircularArc failed");
+            }
+        }
+    }
+
+    // TODO: Select periodic verticies
+    // auto sk.select_periodic_vertex(v0,v1,90.0);
+
+    // Create Mesh
+    Mesh me{sk};
+    me.create();
+
+    me.MaximumElementSize = 0.1;
+    me.MinimumElementSize = 0.01;
+    me.MinimumElementQuality = 0.5;
+
+    result = me.refine();
+    ASSERT_TRUE(result);
+
+    me.save_as(SAVE_DIR, std::string("quarter_circle_mirror_copy_mesh"));
 }
\ No newline at end of file