Implement global objective function data logging for PSO and MPI/PSO

apps/Switched Winding Synchrel/designspace.cpp

@@ -34,7 +34,7 @@ void limit_angular_thickness(CoordinateSpace &d, Particle &p) {
 
 void limit_radial_thickness(CoordinateSpace &d, Particle &p) {
     double_t sum_rt{0.0};
-    for (size_t i = 0; i != 7; ++i) {
+    for (size_t i = 0; i != 6; ++i) {
         std::string key = "rt" + std::to_string(i);
         sum_rt += p[key].Position;
     }
@@ -43,7 +43,7 @@ void limit_radial_thickness(CoordinateSpace &d, Particle &p) {
         for (size_t iter = 0; iter != 3; ++iter) {
             double_t scale{1.0 / sum_rt};
             sum_rt = 0;
-            for (size_t i = 0; i != 7; ++i) {
+            for (size_t i = 0; i != 6; ++i) {
                 std::string key = "rt" + std::to_string(i);
 
                 if (iter == 0) {
@@ -69,6 +69,8 @@ CoordinateSpace design_space() {
 
     double_t inch{25.4e-3};
 
+    d["Jmax"] = Coordinate{30e6,30e6,30e6};
+
     d["ls"] = Coordinate{1*inch,3.5*inch,6*inch};
     d["rro"] = Coordinate{2*inch,3.25*inch,4.5*inch};
 
@@ -77,7 +79,7 @@ CoordinateSpace design_space() {
         d[key] = Coordinate{0.07,0.14,0.21};
     }
 
-    for (size_t i = 0; i != 7; ++i) {
+    for (size_t i = 0; i != 6; ++i) {
         std::string key = "rt" + std::to_string(i);
         d[key] = Coordinate{0.07,0.14,0.21};
     }

apps/Switched Winding Synchrel/globalobjective.cpp

@@ -16,7 +16,7 @@ ObjectiveMap evaluate(Particle particle) {
         omp_set_num_threads(omp_get_max_threads());
     #endif
 
-    size_t Na{5}, Nj{5};
+    size_t Na{11}, Nj{10};
     Eigen::MatrixXd angle = Eigen::MatrixXd::Zero(Na,Nj);
     Eigen::MatrixXd current = Eigen::MatrixXd::Zero(Na,Nj);
     Eigen::MatrixXd torque = Eigen::MatrixXd::Zero(Na,Nj);
@@ -45,7 +45,8 @@ ObjectiveMap evaluate(Particle particle) {
 
     // Series mode torque-speed curve
     //get_torque_and_flux_linkage(model.msp, model.position, Ns, model.pp, torque, flux_cos, flux_sin, angle, current);
-    model.simulate(torque, flux_cos, flux_sin, angle, current);
+    double_t series_conduction_losses;
+    model.simulate(torque, flux_cos, flux_sin, angle, current, series_conduction_losses);
 
     double_t max_torque = torque.array().maxCoeff();
 
@@ -55,7 +56,8 @@ ObjectiveMap evaluate(Particle particle) {
     current.array() *= 0.5;
 
     //get_torque_and_flux_linkage(model.msp, model.position, Ns, model.pp, torque, flux_cos, flux_sin, angle, current);
-    model.simulate(torque, flux_cos, flux_sin, angle, current);
+    double_t parallel_conduction_losses;
+    model.simulate(torque, flux_cos, flux_sin, angle, current, parallel_conduction_losses);
 
     for (size_t i = 0; i != 3; ++i) {
         flux_cos[i].array() *= 0.5;
@@ -89,11 +91,17 @@ ObjectiveMap evaluate(Particle particle) {
 
     om["total volume"] = total_volume;
 
-    //std::cout << min_power << "," << max_torque << "," << active_volume << "," << total_volume << std::endl;
+    om["losses"] = series_conduction_losses;
+
+    std::cout << om["power"] << "," << om["losses"] << std::endl;
 
     return om;
 }
 
+std::vector<std::string> global_objective_keys() {
+    return {"power","torque","active volume","total volume","losses"};
+}
+
 std::function<ObjectiveMap(Particle)> global_objective() {
     return evaluate;
 }
\ No newline at end of file

apps/Switched Winding Synchrel/globalobjective.h

@@ -5,6 +5,8 @@
 
 std::function<ObjectiveMap(Particle)> global_objective();
 
+std::vector<std::string> global_objective_keys();
+
 ObjectiveMap evaluate(Particle p);
 
 #endif // OE_APP_SWITCHED_WINDING_SYNCHREL_GLOBALOBJECTIVE_H
\ No newline at end of file

apps/Switched Winding Synchrel/main.cpp

@@ -13,9 +13,9 @@ int main(int argc, char** argv) {
     MPI_Comm_size(MPI_COMM_WORLD,&size);
 
     // Arguments
-    size_t swarm_size{31};
-    double_t objective_tolerance{1e-2 * 22.0};
-    size_t maximum_iterations{20};
+    size_t swarm_size{35};
+    double_t objective_tolerance{1e-2};
+    size_t maximum_iterations{200};
 
     CoordinateSpace ds = design_space();
 
@@ -25,18 +25,21 @@ int main(int argc, char** argv) {
 
     std::function<ObjectiveMap(Particle)> go = global_objective();
 
-    Swarm swarm{ds,go,swarm_size,objective_tolerance,maximum_iterations,std::greater<double_t>()};
+    Swarm swarm{ds,go,swarm_size,objective_tolerance,maximum_iterations,std::greater<double_t>(),global_objective_keys()};
 
     if (rank == 0) {
         swarm.set_output_file(SAVE_DIR,"sws");
     }
 
-    for (auto &p : swarm.Particles) {
-        p.local_objective() = [](ObjectiveMap m) {
+    for (size_t i = 0; i != swarm_size; ++i) {
+        double_t wi = (i+1.0) / (swarm_size + 1.0);
+        swarm.Particles[i].local_objective() = [wi](ObjectiveMap m) {
+            double_t wpd{(wi) / 9.8};
+            double_t wls{(1.0 - wi) / 0.91};
             if (m["power"] < 55.0) {
-                return ((m["power"] - 55.0) / m["active volume"]);
+                return ((m["power"] - 55.0) / m["total volume"]) / 9.8;
             } else {
-                return m["power"] / m["active volume"];
+                return wpd * m["power"] / m["total volume"] + wls * m["power"] / (m["power"] + m["losses"]);
             }
         };
     }
@@ -49,7 +52,6 @@ int main(int argc, char** argv) {
 
     if (rank == 0) {
         std::cout << swarm;
-        std::cerr << "Reminder: Implement analytical bridge thickness model" << std::endl;
     }
 
     MPI_Finalize();

apps/Switched Winding Synchrel/model.cpp

@@ -129,18 +129,23 @@ void Model::create_physics() {
 void Model::create_post_processor() {
     pp.add<AirgapTorque>(msp,"Torque",fem->region(airgap_contour), Params.rro*2.0/3.0 + Params.rsi*1.0/3.0, Params.rro*1.0/3.0 + Params.rsi*2.0/3.0, 1.0, Params.Np * Params.ls);
     pp.add<Integral>(msp,"FluxA",fem->region({slot_contours[0],slot_contours[1]}),Params.Np * Params.SlotFill * Params.Jmax / Params.Imax * Params.ls, false);
-    pp.add<Integral>(msp,"FluxB",fem->region({slot_contours[2],slot_contours[3]}),Params.Np * Params.SlotFill * Params.Jmax / Params.Imax* Params.ls, false);
-    pp.add<Integral>(msp,"FluxC",fem->region({slot_contours[4],slot_contours[5]}),Params.Np * Params.SlotFill * Params.Jmax / Params.Imax* Params.ls, false);
+    pp.add<Integral>(msp,"FluxB",fem->region({slot_contours[2],slot_contours[3]}),Params.Np * Params.SlotFill * Params.Jmax / Params.Imax * Params.ls, false);
+    pp.add<Integral>(msp,"FluxC",fem->region({slot_contours[4],slot_contours[5]}),Params.Np * Params.SlotFill * Params.Jmax / Params.Imax * Params.ls, false);
+
+    double_t winding_length = Params.ls + M_PI * dws.outer_radius() * sin(M_PI / 2 / Params.Np); // including end turn length estimate
+    double_t resistivity = 1.68e-8 * (1 + 0.00404*(150 - 20));
+    double_t loss_density = resistivity * Params.Jmax * Params.Jmax * Params.SlotFill;
+    pp.add<Integral>(msp,"ConductionLosses",fem->region({slot_contours[0]}), loss_density * Params.Nt * winding_length * 1e-3, false);
 }
 
-void Model::simulate(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>  &flux_cos, std::array<Eigen::MatrixXd,3>  &flux_sin, Eigen::MatrixXd angle, Eigen::MatrixXd current) {
+void Model::simulate(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>  &flux_cos, std::array<Eigen::MatrixXd,3>  &flux_sin, Eigen::MatrixXd angle, Eigen::MatrixXd current, double_t &conduction_losses) {
     torque.setZero();
-    for (size_t i = 0; i!=3; ++i) {
+    for (size_t i = 0; i != 3; ++i) {
         flux_cos[i].setZero();
         flux_sin[i].setZero();
     }
 
-    size_t Ns{16};
+    size_t Ns{32};
     double_t t{0.0};
     double_t dt = 2.0 / (4.0 * position->size());
     size_t dp{position->size() / Ns};
@@ -181,4 +186,6 @@ void Model::simulate(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>  &fl
         *position += dp;
     }
     *position += Ns * dp;
+
+    conduction_losses = pp("ConductionLosses",Oe::VectorXd::Ones(previous[0].size()));
 }

apps/Switched Winding Synchrel/model.h

@@ -49,7 +49,7 @@ public:
     // Post Processor
     PostProcessorInterface pp;
 
-    void simulate(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>  &flux_cos, std::array<Eigen::MatrixXd,3>  &flux_sin, Eigen::MatrixXd angle, Eigen::MatrixXd current);
+    void simulate(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>  &flux_cos, std::array<Eigen::MatrixXd,3>  &flux_sin, Eigen::MatrixXd angle, Eigen::MatrixXd current, double_t &conduction_losses);
 
 protected:
     void create_sketch();

apps/Switched Winding Synchrel/parameters.cpp

@@ -3,14 +3,12 @@
 Parameters::Parameters(Particle p) {
     std::cerr << "Reminder: Implement analytical bridge thickness model" << std::endl;
 
+    Jmax = p["Jmax"].Position;
+
     ls = p["ls"].Position;
     rro = p["rro"].Position;
 
     rsi = rro + 0.03*INCH; // airgap length 30mils
-
-    angular_bridge_thickness = {0.5e-3,0.5e-3,0.5e-3};
-    radial_bridge_thickness = {0.5e-3,0.5e-3,0.5e-3};
-
     std::vector<double_t> rt(6,0.0);
     for (size_t i = 0; i != 6; ++i) {
         std::string key = "rt" + std::to_string(i);
@@ -25,9 +23,130 @@ Parameters::Parameters(Particle p) {
     }
     angular_thickness = at;
 
+    calculate_bridge_thickness();
+
     slot_width = p["sw"].Position * (2 * M_PI * rro) / (Nt);
 
     slot_depth = (p["sd"].Position / p["sw"].Position) * (Imax * Nt) / (2 * M_PI * rro * Jmax * SlotFill);
 
     backiron_thickness = p["bi"].Position * (2 * M_PI * (rro + slot_depth - rro * (1 - p["sw"].Position))) / (2 * Np);
+}
+
+void Parameters::calculate_bridge_thickness() {
+    angular_bridge_thickness = {2*SteelThickness,2*SteelThickness,2*SteelThickness};
+    radial_bridge_thickness = {0.0, 0.0, 0.0};
+
+    auto &rt = radial_thickness;
+    auto &at = angular_thickness;
+
+    for (size_t i = 0; i != 3; ++i) {
+        size_t j = 2 * i;
+        double_t x00 = rt[j] * cos(M_PI / Np);
+        double_t y00 = rt[j] * sin(M_PI / Np);
+        double_t x01 = rro * cos(at[j] * M_PI / 180.0);
+        double_t y01 = rro * sin(at[j] * M_PI / 180.0);
+
+        double_t x10 = rt[j + 1] * cos(M_PI / Np);
+        double_t y10 = rt[j + 1] * sin(M_PI / Np);
+        double_t x11 = rro * cos(at[j + 1] * M_PI / 180.0);
+        double_t y11 = rro * sin(at[j + 1] * M_PI / 180.0);
+
+        Circle c0 = Circle::calculate(x00, y00, x01, y01, Ray{0.0, 0.0, M_PI / Np});
+        Circle c1 = Circle::calculate(x10, y10, x11, y11, Ray{0.0, 0.0, M_PI / Np});
+
+        double_t x0 = c0.x;
+        double_t y0 = c0.y;
+        double_t r0 = c0.r;
+        double_t a00 = atan2(y00 - c0.y, x00 - c0.x);
+        double_t a01 = atan2(y01 - c0.y, x01 - c0.x);
+
+        double_t x1 = c1.x;
+        double_t y1 = c1.y;
+        double_t r1 = c1.r;
+        double_t a10 = atan2(c1.y - y10, c1.x - x10);
+        double_t a11 = atan2(c1.y - y11, c1.x - x11);
+
+        //quadrature rule
+        std::vector<double_t> qp{-1.0 / 3.0 * sqrt(5.0 + 2.0 * sqrt(10.0 / 7.0)),
+                                 -1.0 / 3.0 * sqrt(5.0 - 2.0 * sqrt(10.0 / 7.0)), 0.0,
+                                 +1.0 / 3.0 * sqrt(5.0 - 2.0 * sqrt(10.0 / 7.0)),
+                                 1.0 / 3.0 * sqrt(5.0 + 2.0 * sqrt(10.0 / 7.0))};
+
+        std::vector<double_t> qw{(322 - 13 * sqrt(70)) / 900, (322 + 13 * sqrt(70)) / 900, 128.0 / 225.0,
+                                 (322 + 13 * sqrt(70)) / 900, (322 - 13 * sqrt(70)) / 900};
+
+        for (size_t i = 0; i != qp.size(); ++i) {
+            qp[i] += 1.0;
+            qp[i] /= 2.0;
+            qw[i] /= 2.0;
+        }
+
+        double_t Fx{0.0}, Fy{0.0};
+        for (size_t i = 0; i != qp.size(); ++i) {
+            for (size_t j = 0; j != qp.size(); ++j) {
+                double_t s = qp[i];
+                double_t t = qp[j];
+                double_t w = qw[i] * qw[j];
+
+                double_t xc = x0 * (1 - s) + x1 * s;
+                double_t yc = y0 * (1 - s) + y1 * s;
+                double_t rc = r0 * (1 - s) + r1 * s;
+                double_t a0 = a00 * (1 - s) + a01 * s;
+                double_t a1 = a10 * (1 - s) + a11 * s;
+                double_t ac = a0 * (1 - t) + a1 * t;
+
+                double_t dxcds = x1 - x0;
+                double_t dycds = y1 - y0;
+                double_t drcds = r1 - r0;
+                double_t da0ds = a01 - a00;
+                double_t da1ds = a11 - a10;
+                double_t dacds = da0ds * (1 - t) + da1ds * t;
+                double_t dacdt = a1 - a0;
+
+                double_t x = xc + rc * cos(ac);
+                double_t y = yc + rc * sin(ac);
+
+                double_t dxds = dxcds + drcds * cos(ac) + rc * (-sin(ac) * dacds);
+                double_t dxdt = rc * (-sin(ac) * dacdt);
+                double_t dyds = dycds + drcds * sin(ac) + rc * (+cos(ac) * dacds);
+                double_t dydt = rc * (+cos(ac) * dacdt);
+
+                double_t r = sqrt(x * x + y * y);
+                double_t a = atan2(y, x);
+
+                double_t drdx = +x / r;
+                double_t drdy = +y / r;
+                double_t dadx = +y / (r * r);
+                double_t dady = -x / (r * r);
+
+                double_t drds = drdx * dxds + drdy * dyds;
+                double_t drdt = drdx * dxdt + drdy * dydt;
+                double_t dads = dadx * dxds + dady * dyds;
+                double_t dadt = dadx * dxdt + dady * dydt;
+
+                double_t J = std::abs(drds * dadt - drdt * dads);
+
+                Fx += x * r * J * w;
+                Fy += y * r * J * w;
+            }
+        }
+        Fx *= MaxSpeed * MaxSpeed * SteelDensity;
+        Fy *= MaxSpeed * MaxSpeed * SteelDensity;
+        double_t af = atan(Fy / Fx);
+        double_t Ff = sqrt(Fx * Fx + Fy * Fy);
+        double_t Fr = Ff * cos(M_PI / Np - af);
+
+        for (size_t j = i; j != 3; ++j) {
+            radial_bridge_thickness[j] += Fr / MaximumStress;
+        }
+        radial_bridge_thickness[i] = std::max(radial_bridge_thickness[i], 1.5 * SteelThickness);
+    }
+
+    //radial_bridge_thickness = {2*SteelThickness,2*SteelThickness,2*SteelThickness};
+    // TODO: If mesh is failing, bridge size may be causing errors when it is small
+    /*
+    for (auto v : radial_bridge_thickness) {
+        std::cout << v << std::endl;
+    }
+     */
 }
\ No newline at end of file

apps/Switched Winding Synchrel/parameters.h

@@ -2,6 +2,7 @@
 #define OE_APP_SWITCHED_WINDING_SYNCHREL_PARAMETERS_H
 
 #include "Optimization.hpp"
+#include "Geometry.hpp"
 
 class Parameters {
 public:
@@ -9,8 +10,11 @@ public:
 
     double_t Vdc{600.0}; // dc link voltage
     double_t Imax{400.0*325.0/600.0}; // maximum phase current
-    double_t Jmax{30.0e6}; // maximum current density
+    double_t Jmax{30e6}; // maximum current density
     double_t SlotFill{0.5}; // slot fill factor
+    double_t MaximumStress{200e6};
+    double_t SteelDensity{7.85e3};
+    double_t SteelThickness{0.35e-3};
 
     double_t MaxSpeed{14.0e3*M_PI/30.0}; // rad/s
     double_t CornerSpeed{2.8e3*M_PI/30.0}; // rad/s
@@ -22,20 +26,24 @@ public:
     double_t rri{25.4e-3};
     double_t rro;
     double_t rsi;
+    double_t rso;
 
     std::vector<double_t> angular_bridge_thickness;
     std::vector<double_t> radial_bridge_thickness;
     std::vector<double_t> radial_thickness;
     std::vector<double_t> angular_thickness;
 
-    double_t tooth_face_thickness{2*0.35e-3};
+    double_t tooth_face_thickness{2*SteelThickness};
     double_t slot_depth;
     double_t backiron_thickness;
 
     double_t slot_width;
     double_t slot_opening{1.74e-3};
+
 protected:
     double_t INCH{25.4e-3};
+
+    void calculate_bridge_thickness();
 };
 
 #endif //OERSTED_PARAMETERS_H

src/Optimization/src/CoordinateSpace.cpp

@@ -58,7 +58,7 @@ std::ostream& operator<<(std::ostream& os, CoordinateSpace& p) {
     }
 }
 
-Particle CoordinateSpace::new_particle(std::mt19937 &rng, double_t iv) {
+Particle CoordinateSpace::new_particle(std::mt19937 &rng, double_t iv, std::vector<std::string> go_keys = {}) {
     std::uniform_real_distribution<> dist(-1.0,1.0);
 
     Particle particle;
@@ -74,5 +74,9 @@ Particle CoordinateSpace::new_particle(std::mt19937 &rng, double_t iv) {
         particle.GlobalBest = iv;
     }
 
+    for (auto k : go_keys) {
+        particle.GlobalObjective[k] = NAN;
+    }
+
     return particle;
 }
\ No newline at end of file

src/Optimization/src/CoordinateSpace.h

@@ -11,7 +11,7 @@ class CoordinateSpace {
 public:
     Coordinate &operator[] (std::string key) { return Coordinates[key]; };
 
-    Particle new_particle(std::mt19937 &rng, double_t iv);
+    Particle new_particle(std::mt19937 &rng, double_t iv, std::vector<std::string> go_keys);
 
     size_t size() const { return Coordinates.size(); };
 

src/Optimization/src/Particle.h

@@ -63,6 +63,8 @@ public:
 
     std::function<double_t(ObjectiveMap)> LocalObjective;
 
+    ObjectiveMap GlobalObjective;
+
     double_t PersonalBest{DBL_MAX};
 
     double_t GlobalBest{DBL_MAX};

src/Optimization/src/Swarm.cpp

@@ -7,11 +7,10 @@ void Swarm::run() {
     bool all_converged{false};
     for (size_t iter = 0; iter != MaximumIterations && !all_converged; ++iter) {
         for (size_t i = 0; i != Particles.size(); ++i) {
-            auto go = GlobalObjective(Particles[i]);
-
-            Particles[i].update_personal_best(go,Comparator);
+            Particles[i].GlobalObjective = GlobalObjective(Particles[i]);
+            Particles[i].update_personal_best(Particles[i].GlobalObjective,Comparator);
             for (size_t j = 0; j != Particles.size(); ++j) {
-                Particles[j].update_global_best(go,Particles[i],Comparator);
+                Particles[j].update_global_best(Particles[i].GlobalObjective,Particles[i],Comparator);
             }
         }
 
@@ -57,9 +56,15 @@ void Swarm::set_output_file(std::string path, std::string file_name) {
     fs.open(OutputFile, std::fstream::out);
 
     fs << "iteration,particle,pbest,gbest";
+
+    for (auto key_ind : Particles[0].GlobalObjective.map()) {
+        fs << "," << key_ind.first;
+    }
+
     for (auto key_ind : Particles[0].state().map()) {
         fs << "," << key_ind.first;
     }
+
     fs << "\n";
 
     fs.close();
@@ -71,9 +76,15 @@ void Swarm::append_particle_to_output_file(double_t iter, size_t pnum, Particle
         fs.open(OutputFile, std::fstream::out | std::fstream::app);
 
         fs << iter << "," << pnum << "," << p.personal_best() << "," << p.global_best();
-        for (auto state : p.state().vector()) {
-            fs << "," << state.Position;
+
+        for (auto key_ind : p.GlobalObjective.map()) {
+            fs << "," << p.GlobalObjective[key_ind.first];
         }
+
+        for (auto key_ind : p.state().map()) {
+            fs << "," << p[key_ind.first].Position;
+        }
+
         fs << "\n";
 
         fs.close();
@@ -141,18 +152,19 @@ void Swarm::mpi_master_receive(std::deque<int> &process_queue, std::deque<int> &
         MPI_Iprobe(process, TAG_PARTICLE_INDEX, MPI_COMM_WORLD, &MPI_FLAG, &MPI_STATUS);
         if (MPI_FLAG) {
             int particle;
-            std::vector<double_t> objective(Particles.size(), 0.0);
+            std::vector<double_t> local_objective_values(Particles.size(), 0.0);
 
             MPI_Recv(&particle, 1, MPI_INT, process, TAG_PARTICLE_INDEX, MPI_COMM_WORLD, &MPI_STATUS);
-            MPI_Recv(objective.data(), (int) objective.size(), MPI_DOUBLE, process, TAG_OBJECTIVE_DATA, MPI_COMM_WORLD, &MPI_STATUS);
+            MPI_Recv(local_objective_values.data(), (int) local_objective_values.size(), MPI_DOUBLE, process, TAG_OBJECTIVE_DATA, MPI_COMM_WORLD, &MPI_STATUS);
+            MPI_Recv(Particles[particle].GlobalObjective.vector().data(), (int)Particles[particle].GlobalObjective.vector().size(), MPI_DOUBLE, process, TAG_OBJECTIVE_DATA, MPI_COMM_WORLD, &MPI_STATUS);
 
             // Update global bests
             for (size_t j = 0; j != Particles.size(); ++j) {
-                Particles[j].update_global_best(objective[j],Particles[particle],Comparator);
+                Particles[j].update_global_best(local_objective_values[j],Particles[particle],Comparator);
             }
 
             // Update personal best
-            Particles[particle].update_personal_best(objective[particle],Comparator);
+            Particles[particle].update_personal_best(local_objective_values[particle],Comparator);
 
             append_particle_to_output_file(iter,(size_t)particle,Particles[particle]);
 
@@ -202,6 +214,7 @@ void Swarm::mpi_slave() {
 
             MPI_Send(&particle,1,MPI_INT,0,TAG_PARTICLE_INDEX,MPI_COMM_WORLD);
             MPI_Send(objective.data(),(int)objective.size(),MPI_DOUBLE,0,TAG_OBJECTIVE_DATA,MPI_COMM_WORLD);
+            MPI_Send(go.vector().data(),(int)go.vector().size(),MPI_DOUBLE,0,TAG_OBJECTIVE_DATA,MPI_COMM_WORLD);
         }
     }
 }

src/Optimization/src/Swarm.h

@@ -8,7 +8,7 @@
 
 class Swarm {
 public:
-    Swarm(CoordinateSpace ss, std::function<ObjectiveMap(Particle)> go, size_t Np, double_t otol, size_t maxit, std::function<bool(double_t,double_t)> comp = std::less<double_t>()) : Space{ss}, GlobalObjective{go}, ObjectiveTolerance{otol}, MaximumIterations{maxit}, Comparator{comp} {
+    Swarm(CoordinateSpace ss, std::function<ObjectiveMap(Particle)> go, size_t Np, double_t otol, size_t maxit, std::function<bool(double_t,double_t)> comp = std::less<double_t>(), std::vector<std::string> go_keys = {}) : Space{ss}, GlobalObjective{go}, ObjectiveTolerance{otol}, MaximumIterations{maxit}, Comparator{comp} {
         std::random_device rd;
         RNG = std::mt19937{rd()};
 
@@ -20,7 +20,7 @@ public:
         }
 
         for (size_t i = 0; i != Np; ++i) {
-            Particles.push_back(Space.new_particle(RNG,piv));
+            Particles.push_back(Space.new_particle(RNG,piv,go_keys));
         }
 
         MPI_Initialized(&MPI_FLAG);

test/Model Templates/test_ModelTemplates.hpp

@@ -6,6 +6,7 @@
 #include "ModelTemplates.hpp"
 #include "Physics.hpp"
 #include "Optimization.hpp"
+#include "Geometry.hpp"
 
 #include "gtest.h"
 

test/Model Templates/test_SynchronousReluctanceMachine.cpp

@@ -146,8 +146,7 @@ Parameters calculate_parameters(std::map<std::string,double_t> dsp) {
 
     p.rsi = p.rro + 0.03*25.4e-3; // airgap length 30mils
 
-    p.angular_bridge_thickness = {0.5e-3,0.5e-3,0.5e-3};
-    p.radial_bridge_thickness = {0.5e-3,0.5e-3,0.5e-3};
+    p.angular_bridge_thickness = {0.7e-3,0.7e-3,0.7e-3};
 
     std::vector<double_t> rt(6,0.0);
     for (size_t i = 0; i != 6; ++i) {
@@ -163,6 +162,106 @@ Parameters calculate_parameters(std::map<std::string,double_t> dsp) {
     }
     p.angular_thickness = at;
 
+    p.radial_bridge_thickness = {0.0,0.0,0.0};
+    for (size_t i = 0; i != 3; ++i) {
+        size_t j = 2 * i;
+        double_t x00 = rt[j] * cos(M_PI / p.Np);
+        double_t y00 = rt[j] * sin(M_PI / p.Np);
+        double_t x01 = p.rro * cos(at[j] * M_PI / 180.0);
+        double_t y01 = p.rro * sin(at[j]* M_PI / 180.0);
+
+        double_t x10 = rt[j+1] * cos(M_PI / p.Np);
+        double_t y10 = rt[j+1] * sin(M_PI / p.Np);
+        double_t x11 = p.rro * cos(at[j+1] * M_PI / 180.0);
+        double_t y11 = p.rro * sin(at[j+1] * M_PI / 180.0);
+
+        Circle c0 = Circle::calculate(x00,y00,x01,y01,Ray{0.0,0.0,M_PI / p.Np});
+        Circle c1 = Circle::calculate(x10,y10,x11,y11,Ray{0.0,0.0,M_PI / p.Np});
+
+        double_t x0 = c0.x;
+        double_t y0 = c0.y;
+        double_t r0 = c0.r;
+        double_t a00 = atan2(y00 - c0.y, x00 - c0.x);
+        double_t a01 = atan2(y01 - c0.y, x01 - c0.x);
+
+        double_t x1 = c1.x;
+        double_t y1 = c1.y;
+        double_t r1 = c1.r;
+        double_t a10 = atan2(c1.y - y10, c1.x - x10);
+        double_t a11 = atan2(c1.y - y11, c1.x - x11);
+
+        //quadrature rule
+        std::vector<double_t> qp{-1.0/3.0*sqrt(5.0+2.0*sqrt(10.0/7.0)),-1.0/3.0*sqrt(5.0-2.0*sqrt(10.0/7.0)),0.0,+1.0/3.0*sqrt(5.0-2.0*sqrt(10.0/7.0)),1.0/3.0*sqrt(5.0+2.0*sqrt(10.0/7.0))};
+        std::vector<double_t> qw{(322-13*sqrt(70))/900,(322+13*sqrt(70))/900,128.0/225.0,(322+13*sqrt(70))/900,(322-13*sqrt(70))/900};
+        for (size_t i = 0; i != qp.size(); ++i) {
+            qp[i] += 1.0;
+            qp[i] /= 2.0;