Initial implementation of CylindricalAirgap class which automates creation of...

Initial implementation of CylindricalAirgap class which automates creation of multiple annuli between two radii

src/Sketch/include/Sketch.hpp

@@ -39,4 +39,6 @@
 
 #include "../src/Contour.h"
 
+#include "../src/CylindricalAirgap.h"
+
 #endif
\ No newline at end of file

src/Sketch/src/CylindricalAirgap.cpp

+#include "CylindricalAirgap.h"
\ No newline at end of file

src/Sketch/src/CylindricalAirgap.h

+#ifndef OERSTED_CYLINDRICALAIRGAP_H
+#define OERSTED_CYLINDRICALAIRGAP_H
+
+#include <cmath>
+#include "Sketch.h"
+#include "LineSegment.h"
+#include "CircularArc.h"
+#include "Coincident.h"
+#include "Angle.h"
+
+// TODO: Move to a library for higher order functions
+
+template <size_t Nl>
+class CylindricalAirgap {
+public:
+    CylindricalAirgap() {};
+
+    CylindricalAirgap(Sketch &sketch, std::shared_ptr<const Vertex> vc, std::shared_ptr<const Vertex> vi0, std::shared_ptr<const Vertex> vo0, double_t angle) {
+        sketch.solve();
+
+        // Get matching periodic verticies
+        auto vi1 = sketch.select_periodic_vertex(vi0, vc, angle);
+        auto vo1 = sketch.select_periodic_vertex(vo0, vc, angle);
+
+        // Determine radii
+        double_t ri = vc->hypot(vi0);
+        double_t ro = vc->hypot(vo0);
+
+        double_t ai = vc->atan(vi0);
+        double_t ao = vc->atan(vo0);
+
+        if (abs(ai - ao) > FLT_EPSILON * M_PI) { // TODO: Check for branch cut discontinuity due to numerical errors
+            std::cerr << "Periodic boundaries should be aligned when creating airgap. Airgap will not be created" << std::endl;
+            return;
+        }
+
+        double_t aa = (ai + ao) / 2.0;
+        double_t da = angle * M_PI / 180.0;
+
+        double_t r = ri;
+        double_t dr = (ro - ri) / (Nl + 1);
+
+        auto v0 = vi0;
+        auto v1 = vi1;
+        for (size_t i = 0; i != Nl; ++i) {
+            r += dr;
+
+            // Create Airgap Vertices
+            Vertex0[i] = sketch.new_element<Vertex>(r * std::cos(aa), r * std::sin(aa));
+            Vertex1[i] = sketch.new_element<Vertex>(r * std::cos(aa + da), r * std::sin(aa + da));
+
+            // Create Airgap Edge LineSegments
+            Line0[i] = sketch.new_element<LineSegment>(v0, Vertex0[i]);
+            Line1[i] = sketch.new_element<LineSegment>(v1, Vertex1[i]);
+
+            // Create Airgap CircularArc Interface
+            Interface[i] = sketch.new_element<CircularArc>(Vertex0[i], Vertex1[i], vc, r);
+
+            // Next
+            v0 = Vertex0[i];
+            v1 = Vertex1[i];
+        }
+
+        // Last line segments
+        Line0[Nl+1] = sketch.new_element<LineSegment>(v0,vo0);
+        Line1[Nl+1] = sketch.new_element<LineSegment>(v1,vo1);
+
+            /*
+            // Create Airgap Vertices
+            Vertex0 = sketch.new_element<Vertex>(ra * std::cos(aa), ra * std::sin(aa));
+            Vertex1 = sketch.new_element<Vertex>(ra * std::cos(aa + da), ra*std::sin(aa + da));
+
+            // Create Airgap Edge LineSegments
+            Line0ia = sketch.new_element<LineSegment>(vi0, Vertex0);
+            Line0ao = sketch.new_element<LineSegment>(Vertex0, vo0);
+
+            Line1ia = sketch.new_element<LineSegment>(vi1, Vertex1);
+            Line1ao = sketch.new_element<LineSegment>(Vertex1, vo1);
+
+            // Create Airgap CircularArc Interface
+            AirgapInterface = sketch.new_element<CircularArc>(Vertex0, Vertex1, vc, ra);
+
+            // TODO: Should constraints be added?
+             */
+    }
+
+    std::shared_ptr<Vertex> vertex0(size_t i) const { return Vertex0[i]; };
+
+    std::shared_ptr<Vertex> vertex1(size_t i) const { return Vertex1[i]; };
+
+    std::shared_ptr<LineSegment> line0(size_t i) const { return Line0[i]; };
+
+    std::shared_ptr<LineSegment> line1(size_t i) const { return Line1[i]; };
+
+    std::shared_ptr<CircularArc> interface(size_t i) const { return Interface[i]; };
+
+protected:
+    std::array<std::shared_ptr<Vertex>,Nl> Vertex0;
+    std::array<std::shared_ptr<Vertex>,Nl> Vertex1;
+
+    std::array<std::shared_ptr<LineSegment>,Nl+1> Line0;
+    std::array<std::shared_ptr<LineSegment>,Nl+1> Line1;
+
+    std::array<std::shared_ptr<CircularArc>,Nl> Interface;
+};
+
+#endif //OERSTED_CYLINDRICALAIRGAP_H

src/Sketch/src/Vertex.cpp

@@ -6,7 +6,7 @@ void Vertex::register_parameters(Sketch *s) const {
     s->add_parameter(std::const_pointer_cast<Variable>(Y));
 };
 
-double2 Vertex::rotate(std::shared_ptr<Vertex const> const &origin, double angle) const {
+double2 Vertex::rotate(std::shared_ptr<Vertex const> origin, double angle) const {
     double x = origin->x();
     double y = origin->y();
 

src/Sketch/src/Vertex.h

@@ -9,11 +9,11 @@ class Vertex : public SketchElement {
 public:
     Vertex() : X(std::make_shared<Variable const>()), Y(std::make_shared<Variable const>()) {};
 
-    Vertex(double x, double y) : X(std::make_shared<Variable const>(x)), Y(std::make_shared<Variable const>(y)) {};
+    Vertex(double_t x, double_t y) : X(std::make_shared<Variable const>(x)), Y(std::make_shared<Variable const>(y)) {};
 
     Vertex(std::shared_ptr<Variable const> x, std::shared_ptr<Variable const> y) : X(x), Y(y) {};
 
-    Vertex(std::shared_ptr<Vertex const> const &v) : X(v->X), Y(v->Y) {};
+    Vertex(std::shared_ptr<Vertex const> v) : X(v->X), Y(v->Y) {};
 
     size_t set_equation_index(size_t i) override {
         EquationIndex = i;
@@ -26,7 +26,7 @@ public:
 
     bool operator==(Vertex const &v) { return (v.X == X) && (v.Y == Y); };
 
-    bool is_identical(std::shared_ptr<Vertex const> const &v, std::shared_ptr<Vertex const> const &origin, double_t angle) const {
+    bool is_identical(std::shared_ptr<Vertex const> v, std::shared_ptr<Vertex const> origin, double_t angle) const {
         double_t xr = v->x() - origin->x();
         double_t yr = v->y() - origin->y();
         double_t rr = sqrt(xr * xr + yr * yr);
@@ -39,19 +39,23 @@ public:
         return (sqrt(xr * xr + yr * yr) < rr * FLT_EPSILON);
     }
 
-    double atan() const { return std::atan2(y(), x()); };
+    double_t atan() const { return std::atan2(y(), x()); };
 
-    double hypot() const { return std::hypot(x(), y()); };
+    double_t hypot() const { return std::hypot(x(), y()); };
 
-    double x() const { return X->value(); };
+    double_t hypot(std::shared_ptr<const Vertex> v) const { return std::hypot(x() - v->x(), y() - v->y()); };
 
-    double y() const { return Y->value(); };
+    double_t atan(std::shared_ptr<const Vertex> v) const { return std::atan2(v->y() - y(), v->x() - x()); };
+
+    double_t x() const { return X->value(); };
+
+    double_t y() const { return Y->value(); };
 
     void register_parameters(Sketch *s) const override;
 
     void update(Eigen::MatrixXd &J, Eigen::VectorXd &r) const override {};
 
-    double2 rotate(std::shared_ptr<Vertex const> const &origin, double angle) const;
+    double2 rotate(std::shared_ptr<const Vertex> origin, double_t angle) const;
 
 protected:
     std::shared_ptr<Variable const> X;

test/LibraryIntegration/Mesh_To_FEM_Test.cpp

@@ -913,6 +913,7 @@ public:
         auto mr0 = sk.new_element<MirrorCopy>(mirror_vec_r, lr3);
 
         // Rotor part of airgap
+        /*
         auto vra0 = sk.new_element<Vertex>(rro + airgap / 2.0, 0.0);
         auto vra1 = sk.new_element<Vertex>(0.0, rro + airgap / 2.0);
 
@@ -924,6 +925,7 @@ public:
         auto lra1 = sk.new_element<LineSegment>(vr1p, vra1);
 
         continuous_airgap = sk.new_element<CircularArc>(vra0, vra1, origin, rro + airgap / 2.0);
+        */
 
         // Stator
         rsi = rro + airgap;
@@ -979,12 +981,31 @@ public:
         residual = sk.solve();
         EXPECT_LE(residual, FLT_EPSILON * rro);
 
+        // Rotor part of airgap
+        std::cout << vr1->x() << "," << vr1->y() << std::endl;
+        std::cout << vs0->x() << "," << vs0->y() << std::endl;
+
+        CylindricalAirgap<1> cyag{sk,origin,vr1,vs0,90.0};
+
+        //auto vra0 = sk.new_element<Vertex>(rro + airgap / 2.0, 0.0);
+        //auto vra1 = sk.new_element<Vertex>(0.0, rro + airgap / 2.0);
+
+        auto fra0 = sk.new_element<Fixation>(cyag.vertex0(0));
+        auto fra1 = sk.new_element<Fixation>(cyag.vertex1(0));
+
+        //auto vr1p = sk.select_periodic_vertex(vr1, origin, 90.0);
+        //auto lra0 = sk.new_element<LineSegment>(vr1, vra0);
+        //auto lra1 = sk.new_element<LineSegment>(vr1p, vra1);
+
+        //continuous_airgap = sk.new_element<CircularArc>(vra0, vra1, origin, rro + airgap / 2.0);
+        continuous_airgap = cyag.interface(0);
+
         // Stator part of airgap
-        auto vsa0 = vs0;
-        auto vsa1 = sk.select_periodic_vertex(vsa0, origin, 90.0);
+        //auto vsa0 = vs0;
+        //auto vsa1 = sk.select_periodic_vertex(vsa0, origin, 90.0);
 
-        auto lsa0 = sk.new_element<LineSegment>(vra0, vsa0);
-        auto lsa1 = sk.new_element<LineSegment>(vra1, vsa1);
+        //auto lsa0 = sk.new_element<LineSegment>(vra0, vsa0);
+        //auto lsa1 = sk.new_element<LineSegment>(vra1, vsa1);
 
         // Build and Save
         residual = sk.solve();