Update to SynchronousReluctanceRotor to improve geometry quality for certain parameter combinations

Improvements to newton raphson convergence checking which improves computation time

Changes IntegralAverage post processor to Integral with an optional switch to calculate the average

Begin test implementations of functions for torque-speed and power-speed mapping of synchronous reluctance example

Being test implementation of functions for particle swarm optimization

src/Model Templates/src/SynchronousReluctanceRotor.cpp

@@ -63,8 +63,6 @@ void SynchronousReluctanceRotor::add_to_sketch(Sketch &s, std::shared_ptr<Vertex
 
         auto l0 = s.new_element<LineSegment>(v0,v1);
         L.push_back(l0);
-        //auto l1 = s.new_element<LineSegment>(origin,v2); // for construction
-        //auto l2 = s.new_element<LineSegment>(origin,v3); // for construction
 
         auto c0 = s.new_element<CircularArc>(v2,v3,origin,ro-dra[j]);
         C.push_back(c0);
@@ -73,26 +71,35 @@ void SynchronousReluctanceRotor::add_to_sketch(Sketch &s, std::shared_ptr<Vertex
             Circle c = Circle::calculate(xr[j],yr[j],xa[j],ya[j],Ray{0.0,0.0,M_PI/Np});
 
             auto vc1 = s.new_element<Vertex>(c.x, c.y);
-            s.new_element<Coincident<LineSegment>>(vc1, L[1]);
-
             auto c1 = s.new_element<CircularArc>(v1, v3, vc1, c.r);
             C.push_back(c1);
+
+            double_t ao3 = atan2(v3->y() - origin->y(), v3->x() - origin->x()) + M_PI_2;
+            double_t a3c = atan2(vc1->y() - v3->y(), vc1->x() - v3->x());
+            if (a3c < ao3) {
+                s.new_element<Coincident<LineSegment>>(vc1, L[1]);
+            } else {
+                auto ltana = s.new_element<LineSegment>(origin,v3,true);
+                s.new_element<Tangency>(c1,ltana);
+            }
         }
 
         {
             Circle c = Circle::calculate(xr[j+1],yr[j+1],xa[j+1],ya[j+1],Ray{0.0,0.0,M_PI/Np});
 
             auto vc1 = s.new_element<Vertex>(c.x,c.y);
-            s.new_element<Coincident<LineSegment>>(vc1, L[1]);
-
             auto c1 = s.new_element<CircularArc>(v0, v2, vc1, c.r);
             C.push_back(c1);
-        }
-        //auto vc2 = s.new_element<Vertex>(1.1 * ro * cos(a180p), 1.1 * ro * sin(a180p));
-        //s.new_element<Coincident<LineSegment>>(vc2,L[1]);
 
-        //auto c2 = s.new_element<CircularArc>(v0,v2,vc2,RadialThickness[j]);
-        //C.push_back(c2);
+            double_t ao2 = atan2(v2->y() - origin->y(), v2->x() - origin->x()) + M_PI_2;
+            double_t a2c = atan2(vc1->y() - v2->y(), vc1->x() - v2->x());
+            if (a2c < ao2) {
+                s.new_element<Coincident<LineSegment>>(vc1, L[1]);
+            } else {
+                auto ltan = s.new_element<LineSegment>(origin,v2,true);
+                s.new_element<Tangency>(c1,ltan);
+            }
+        }
     }
 
     // Impose constraints

src/Physics/src/Physics.h

@@ -64,8 +64,10 @@ public:
 
         size_t iter = 0;
         double_t norm_dv = 1.0;
+        double_t norm_r = 1.0;
+        double_t tol = 1e-2;
 
-        while (iter < 30 && norm_dv > sqrt(DBL_EPSILON) * std::sqrt(s->v.size())) {
+        while (iter < 30 && (norm_dv > tol || norm_r > tol)) {
             s->Solver.compute(s->J);
             if (s->Solver.info() != Oe::Success) {
                 std::cerr << "Factorization of Jacobian failed";
@@ -81,6 +83,8 @@ public:
 
             linearize(s);
 
+            norm_r = s->r.norm() / (s->J * s->v).norm();
+
             iter++;
         }
     }

src/Physics/src/PostProcessorInterface.cpp

@@ -5,7 +5,7 @@ void PostProcessorInterface::register_element(std::shared_ptr<PhysicsInterface>
     pi->register_post_processor_element(ppe);
 }
 
-void IntegralAverage::build() {
+void Integral::build() {
     // Determines the matrix representing the discrete integral average operator of the unknown field over the given regions
     auto pptr = Physics.lock();
 
@@ -19,14 +19,12 @@ void IntegralAverage::build() {
         }
     }
 
-    double_t area{(weights(0) * indicator).sum()}; // Area = integral(1.0*dOmega)
+    Area = (weights(0) * indicator).sum(); // Area = integral(1.0*dOmega)
     OperatorMatrix = basis(0) * (weights(0) * indicator).matrix();
     for (size_t i = 1; i != basis.size(); ++i) {
-        area += (weights(i) * indicator).sum();
+        Area += (weights(i) * indicator).sum();
         OperatorMatrix += basis(i) * (weights(i) * indicator).matrix();
     }
-
-    OperatorMatrix /= area;
 }
 
 double_t LorentzTorque::operator()(Oe::VectorXd const &v) {

src/Physics/src/PostProcessorInterface.h

@@ -22,20 +22,34 @@ protected:
     std::weak_ptr<PhysicsInterface> Physics;
 };
 
-class IntegralAverage : public PostProcessorElement {
+class Integral : public PostProcessorElement {
 public:
-    IntegralAverage(std::shared_ptr<PhysicsInterface> p, std::vector<std::shared_ptr<DiscreteRegion>> r) : PostProcessorElement{p}, Regions{r} {};
+    Integral(std::shared_ptr<PhysicsInterface> p, std::vector<std::shared_ptr<DiscreteRegion>> r, double_t s = 1.0, bool average = false) : PostProcessorElement{p}, Regions{r}, Scale{s}, Average{average} {};
 
-    IntegralAverage(std::shared_ptr<PhysicsInterface> p, std::shared_ptr<DiscreteRegion> r) : PostProcessorElement{p}, Regions{r} {};
+    Integral(std::shared_ptr<PhysicsInterface> p, std::shared_ptr<DiscreteRegion> r, double_t s = 1.0, bool average = false) : PostProcessorElement{p}, Regions{r}, Scale{s}, Average{average} {};
 
     void build() override;
 
-    double_t operator()(Oe::VectorXd const &v) override { return OperatorMatrix.dot(v); };
+    double_t operator()(Oe::VectorXd const &v) override {
+        double_t value = Scale * OperatorMatrix.dot(v);
+
+        if (Average) {
+            value /= Area;
+        }
+
+        return value;
+    };
 
 protected:
     std::vector<std::shared_ptr<DiscreteRegion>> Regions;
 
     Oe::VectorXd OperatorMatrix;
+
+    double_t Scale{1.0};
+
+    bool Average{false};
+
+    double_t Area{0.0};
 };
 
 class LorentzTorque : public PostProcessorElement {
@@ -81,6 +95,10 @@ public:
         register_element(physics, Elements[key]);
     }
 
+    void erase(std::string key) {
+        Elements.erase(key);
+    }
+
     void register_element(std::shared_ptr<PhysicsInterface> pi, std::shared_ptr<PostProcessorElement> ppi);
 
     void build() {
@@ -89,6 +107,8 @@ public:
         }
     }
 
+    // TODO: implement template<T> T get(std::string key) { ... }; where it is attempted to cast Elements[key] to type T and return the shared_ptr if successful, this will allow modification of postprocessor elements
+
 protected:
     std::unordered_map<std::string, std::shared_ptr<PostProcessorElement>> Elements;
 };

test/LibraryIntegration/Mesh_To_FEM_Test.cpp

@@ -1172,7 +1172,7 @@ TEST_F(Salient_Pole_Synchrel, Test) {
         v_prev = solution->v;
 
         //Verify equation is solved
-        EXPECT_LE(solution->r.norm(), FLT_EPSILON * solution->f.norm() * sqrt(solution->f.size()));
+        EXPECT_LE(solution->r.norm(), 1e-2 * (solution->J * solution->v).norm());
 
         // Save image
         Eigen::VectorXd vv = msph->recover_boundary(solution->v);
@@ -1250,7 +1250,7 @@ TEST_F(Salient_Pole_Synchrel, Test_OpenMP) {
             msph.solve(solution, fargs);
 
             //Verify equation is solved
-            EXPECT_LE(solution->r.norm(), FLT_EPSILON * solution->f.norm() * sqrt(solution->f.size()));
+            EXPECT_LE(solution->r.norm(), 1e-2 * (solution->J * solution->v).norm());
         }
 
         // Increment Position

test/Model Templates/test_ModelTemplates.hpp

@@ -8,6 +8,8 @@
 
 #include "gtest.h"
 
+#include <random>
+
 #define SAVE_DIR "./test/output/Model Templates/"
 
 #endif //OERSTED_TEST_MODELTEMPLATES_HPP
\ No newline at end of file

test/Physics/Finite_Element_Mesh_Test.cpp

@@ -681,7 +681,7 @@ TEST_F(TwoRegionHexagon, magnetostatic_forcing_air) {
     // Construct Post-Processor
     PostProcessorInterface pp;
 
-    pp.add<IntegralAverage>(ms, "FluxLinkage", inner_region);
+    pp.add<Integral>(ms, "FluxLinkage", inner_region, 1.0, true);
     pp.build();
 
     Eigen::VectorXd ones = Eigen::VectorXd::Zero(solution->v.size());
@@ -1345,8 +1345,8 @@ TEST_F(Two_Region_Quarter_Octagon_Discontinuous, magnetostatic_periodic) {
 
     // Prepare post-processor
     PostProcessorInterface pp;
-    pp.add<IntegralAverage>(ms, "InnerFluxLinkage", inner_region);
-    pp.add<IntegralAverage>(ms, "OuterFluxLinkage", outer_region);
+    pp.add<Integral>(ms, "InnerFluxLinkage", inner_region, 1.0, true);
+    pp.add<Integral>(ms, "OuterFluxLinkage", outer_region, 1.0, true);
 
     for(size_t iter = 0; iter!=5; ++iter) {
         ms->assemble();
@@ -1416,8 +1416,8 @@ TEST_F(Two_Region_Quarter_Octagon_Discontinuous, magnetostatic_antiperiodic) {
 
     // Prepare post-processor
     PostProcessorInterface pp;
-    pp.add<IntegralAverage>(ms, "InnerFluxLinkage", inner_region);
-    pp.add<IntegralAverage>(ms, "OuterFluxLinkage", outer_region);
+    pp.add<Integral>(ms, "InnerFluxLinkage", inner_region, 1.0, true);
+    pp.add<Integral>(ms, "OuterFluxLinkage", outer_region, 1.0, true);
 
     std::vector<double_t> vv_expected{0.0,0.0,6.60365e-8,0.0,0.0,6.60365e-8,0.0,0.0,0.0,0.0};
     std::vector<double_t> of_expected{2.201e-8,0.0,-2.201e-8,0.0,2.201e-8};