Update Switched Winding Synchrel application to consider losses and torque ripple

Switch to STATIC libraries for building on Eos

CMakeLists.txt

@@ -4,13 +4,14 @@ project(Oersted)
 
 find_package(MPI REQUIRED)
 find_package(OpenMP REQUIRED)
+find_package(Boost REQUIRED)
 
 include_directories(${MPI_INCLUDE_PATH})
 include_directories(./lib/)
 include_directories(./lib/Eigen/)
 include_directories(./lib/GoogleTest/googletest/include/gtest/)
 
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++17")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++14")
 #set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O3 --coverage -fopenmp")
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O3 -fopenmp")
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")

apps/Switched Winding Synchrel/globalobjective.cpp

@@ -2,12 +2,6 @@
 #include "model.h"
 #include "torquespeedcurve.h"
 
-//#include "parameters.h"
-//#include "create_sketch.h"
-//#include "create_mesh.h"
-//#include "create_physics.h"
-//#include "create_post_processor.h"
-
 ObjectiveMap evaluate(Particle particle) {
     Model model{particle};
 
@@ -20,6 +14,7 @@ ObjectiveMap evaluate(Particle particle) {
     Eigen::MatrixXd angle = Eigen::MatrixXd::Zero(Na,Nj);
     Eigen::MatrixXd current = Eigen::MatrixXd::Zero(Na,Nj);
     Eigen::MatrixXd torque = Eigen::MatrixXd::Zero(Na,Nj);
+    Eigen::MatrixXd ripple = Eigen::MatrixXd::Zero(Na,Nj);
     std::array<Eigen::MatrixXd,3> flux_cos;
     std::array<Eigen::MatrixXd,3> flux_sin;
     for (size_t i = 0; i != 3; ++i) {
@@ -27,6 +22,11 @@ ObjectiveMap evaluate(Particle particle) {
         flux_sin[i] = Eigen::MatrixXd::Zero(Na,Nj);
     }
 
+    std::array<Eigen::MatrixXd,4> losses;
+    for (size_t i = 0; i != losses.size(); ++i) {
+        losses[i] = Eigen::MatrixXd::Zero(Na,Nj);
+    }
+
     std::vector<Oe::VectorXd> previous(Na * Nj);
 
     for (size_t i = 0; i != Na; ++i) {
@@ -39,46 +39,41 @@ ObjectiveMap evaluate(Particle particle) {
     // Torque-speed curve values
     Eigen::ArrayXd SPEED = Eigen::ArrayXd::Zero(113);
     Eigen::ArrayXd TORQUE = Eigen::ArrayXd::Zero(113);
+    Eigen::ArrayXd RIPPLE = Eigen::ArrayXd::Zero(113);
+    Eigen::ArrayXd LOSSES = Eigen::ArrayXd::Zero(113);
     for (size_t i = 0; i != SPEED.size(); ++i) {
         SPEED(i) = (model.Params.CornerSpeed + i * (model.Params.MaxSpeed - model.Params.CornerSpeed) / (SPEED.size() - 1)) * 0.030 / M_PI;
     }
 
     // Series mode torque-speed curve
     //get_torque_and_flux_linkage(model.msp, model.position, Ns, model.pp, torque, flux_cos, flux_sin, angle, current);
-    double_t series_conduction_losses;
-    model.simulate(torque, flux_cos, flux_sin, angle, current, series_conduction_losses);
-
+    model.simulate(torque, ripple, flux_cos, flux_sin, angle, current, losses);
     double_t max_torque = torque.array().maxCoeff();
 
-    torque_speed_curve(torque,flux_cos,flux_sin,model.Params.Vdc,model.Params.Np,SPEED,TORQUE);
+    torque_speed_curve(torque,ripple,losses,flux_cos,flux_sin,model.Params.Vdc,model.Params.Np,SPEED,TORQUE,RIPPLE,LOSSES);
 
     // Parallel mode torque-speed curve
     current.array() *= 0.5;
 
     //get_torque_and_flux_linkage(model.msp, model.position, Ns, model.pp, torque, flux_cos, flux_sin, angle, current);
-    double_t parallel_conduction_losses;
-    model.simulate(torque, flux_cos, flux_sin, angle, current, parallel_conduction_losses);
+    model.simulate(torque, ripple, flux_cos, flux_sin, angle, current, losses);
 
     for (size_t i = 0; i != 3; ++i) {
         flux_cos[i].array() *= 0.5;
         flux_sin[i].array() *= 0.5;
     }
 
-    torque_speed_curve(torque,flux_cos,flux_sin,model.Params.Vdc,model.Params.Np,SPEED,TORQUE);
+    torque_speed_curve(torque,ripple,losses,flux_cos,flux_sin,model.Params.Vdc,model.Params.Np,SPEED,TORQUE,RIPPLE,LOSSES);
 
     // Output
     Eigen::ArrayXd POWER = TORQUE * SPEED * M_PI / 30.0;
-    /*
-    for (size_t i = 0; i != SPEED.size(); ++i) {
-        std::cout << SPEED[i] << "," << (size_t)TORQUE[i] << "," << (size_t)POWER[i] << std::endl;
-    }
-    */
 
     // Post Process
     ObjectiveMap om;
 
     om["power"] = POWER.minCoeff();
     om["torque"] = max_torque;
+    om["ripple"] = RIPPLE.maxCoeff();
 
     double_t backiron_radius = model.dws.outer_radius();
     double_t active_volume = model.Params.ls * pow(backiron_radius,2.0) * M_PI * 1e3;
@@ -91,15 +86,15 @@ ObjectiveMap evaluate(Particle particle) {
 
     om["total volume"] = total_volume;
 
-    om["losses"] = series_conduction_losses;
+    om["losses"] = LOSSES.maxCoeff();
 
-    std::cout << om["power"] << "," << om["losses"] << std::endl;
+    std::cout << om["torque"] << "," << om["ripple"] << "," << om["power"] << "," << om["losses"] << std::endl;
 
     return om;
 }
 
 std::vector<std::string> global_objective_keys() {
-    return {"power","torque","active volume","total volume","losses"};
+    return {"power","torque","ripple","active volume","total volume","losses"};
 }
 
 std::function<ObjectiveMap(Particle)> global_objective() {

apps/Switched Winding Synchrel/main.cpp

@@ -13,7 +13,8 @@ int main(int argc, char** argv) {
     MPI_Comm_size(MPI_COMM_WORLD,&size);
 
     // Arguments
-    size_t swarm_size{35};
+    size_t N{11};
+    size_t swarm_size{((N - 2) * (N - 1)) / 2};
     double_t objective_tolerance{1e-2};
     size_t maximum_iterations{200};
 
@@ -31,15 +32,34 @@ int main(int argc, char** argv) {
         swarm.set_output_file(SAVE_DIR,"sws");
     }
 
-    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};
+    double_t dw{1.0/(N-1.0)};
+    double_t wp{0.0},wl{0.0},wr{0.0};
+    for (size_t i = 0; i < swarm_size; ++i) {
+        wp += dw;
+        if (wp + wl > 1.0 - dw + FLT_EPSILON) {
+            wp = dw;
+            wl += dw;
+        }
+        wr = 1.0 - wp - wl;
+
+        if (rank == 0) {
+        std::cout << wp << "," << wl << "," << wr << std::endl;
+        }
+
+        swarm.Particles[i].local_objective() = [wp,wl,wr](ObjectiveMap m) {
+            double_t wpd{wp / 9.8};
+            double_t wls{wl * 5.5};
+            double_t wrl{wr * 0.1};
             if (m["power"] < 55.0) {
-                return ((m["power"] - 55.0) / m["total volume"]) / 9.8;
+                double_t value{0.0};
+                value += (std::min(0.0,m["power"] - 55.0)) / 55.0;
+                return value;
             } else {
-                return wpd * m["power"] / m["total volume"] + wls * m["power"] / (m["power"] + m["losses"]);
+                double_t value{0.0};
+                value += wpd * m["power"] / m["total volume"];
+                value += wls / m["losses"];
+                value += wrl * (m["torque"] / (m["torque"] - m["ripple"]));
+                return value;
             }
         };
     }

apps/Switched Winding Synchrel/model.cpp

@@ -134,23 +134,36 @@ void Model::create_post_processor() {
 
     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);
+    pp.add<Integral>(msp,"ConductionLosses",fem->region({slot_contours[0]}), resistivity * 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, double_t &conduction_losses) {
+void Model::simulate(Eigen::MatrixXd &torque, Eigen::MatrixXd &ripple, std::array<Eigen::MatrixXd,3>  &flux_cos, std::array<Eigen::MatrixXd,3>  &flux_sin, Eigen::MatrixXd angle, Eigen::MatrixXd current, std::array<Eigen::MatrixXd,4> &losses) {
     torque.setZero();
+
     for (size_t i = 0; i != 3; ++i) {
         flux_cos[i].setZero();
         flux_sin[i].setZero();
     }
 
+    for (size_t i = 0; i != 3; ++i) {
+        losses[i].setZero();
+    }
+
+    Eigen::MatrixXd torque_min = +FLT_MAX*Eigen::MatrixXd::Ones(torque.rows(),torque.cols());
+    Eigen::MatrixXd torque_max = -FLT_MAX*Eigen::MatrixXd::Ones(torque.rows(),torque.cols());
+
     size_t Ns{32};
     double_t t{0.0};
     double_t dt = 2.0 / (4.0 * position->size());
     size_t dp{position->size() / Ns};
 
     std::vector<Eigen::VectorXd> previous(angle.size());
+    std::vector<std::array<Eigen::VectorXd,3>> Bx_cos(angle.size());
+    std::vector<std::array<Eigen::VectorXd,3>> By_cos(angle.size());
+
+    std::vector<std::array<Eigen::VectorXd,3>> Bx_sin(angle.size());
+    std::vector<std::array<Eigen::VectorXd,3>> By_sin(angle.size());
+
     for (size_t is = 0; is != Ns; ++is) {
         msp->assemble();
 
@@ -169,9 +182,14 @@ void Model::simulate(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>  &fl
             previous[i] = solution->v;
 
             //Verify equation is solved
-            //EXPECT_LE(solution->r.norm(), 1e-2 * (solution->J * solution->v).norm());
+            double_t tau = pp("Torque", solution->v);
+            torque(i) += tau / Ns;
+            torque_max(i) = std::max(torque_max(i),tau);
+            torque_min(i) = std::min(torque_min(i),tau);
+
+            double_t cosa = (2.0 * cos(2.0 * M_PI * t) / Ns);
+            double_t sina = (2.0 * sin(2.0 * M_PI * t) / Ns);
 
-            torque(i) += pp("Torque", solution->v) / Ns;
             flux_cos[0](i) += pp("FluxA",solution->v) * (2.0 * cos(2.0 * M_PI * t) / Ns);
             flux_cos[1](i) += pp("FluxB",solution->v) * (2.0 * cos(2.0 * M_PI * t) / Ns);
             flux_cos[2](i) += pp("FluxC",solution->v) * (2.0 * cos(2.0 * M_PI * t) / Ns);
@@ -179,6 +197,39 @@ void Model::simulate(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>  &fl
             flux_sin[0](i) += pp("FluxA",solution->v) * (2.0 * sin(2.0 * M_PI * t) / Ns);
             flux_sin[1](i) += pp("FluxB",solution->v) * (2.0 * sin(2.0 * M_PI * t) / Ns);
             flux_sin[2](i) += pp("FluxC",solution->v) * (2.0 * sin(2.0 * M_PI * t) / Ns);
+
+            //
+            auto d = msp->derivatives();
+
+            if (is == 0) {
+                Bx_cos[i][0] = msp->derivatives().dy(0).transpose() * solution->v * cosa;
+                By_cos[i][0] = -msp->derivatives().dx(0).transpose() * solution->v * cosa;
+                Bx_cos[i][1] = msp->derivatives().dy(1).transpose() * solution->v * cosa;
+                By_cos[i][1] = -msp->derivatives().dx(1).transpose() * solution->v * cosa;
+                Bx_cos[i][2] = msp->derivatives().dy(2).transpose() * solution->v * cosa;
+                By_cos[i][2] = -msp->derivatives().dx(2).transpose() * solution->v * cosa;
+
+                Bx_sin[i][0] = msp->derivatives().dy(0).transpose() * solution->v * sina;
+                By_sin[i][0] = -msp->derivatives().dx(0).transpose() * solution->v * sina;
+                Bx_sin[i][1] = msp->derivatives().dy(1).transpose() * solution->v * sina;
+                By_sin[i][1] = -msp->derivatives().dx(1).transpose() * solution->v * sina;
+                Bx_sin[i][2] = msp->derivatives().dy(2).transpose() * solution->v * sina;
+                By_sin[i][2] = -msp->derivatives().dx(2).transpose() * solution->v * sina;
+            } else {
+                Bx_cos[i][0] += msp->derivatives().dy(0).transpose() * solution->v * cosa;
+                By_cos[i][0] += -msp->derivatives().dx(0).transpose() * solution->v * cosa;
+                Bx_cos[i][1] += msp->derivatives().dy(1).transpose() * solution->v * cosa;
+                By_cos[i][1] += -msp->derivatives().dx(1).transpose() * solution->v * cosa;
+                Bx_cos[i][2] += msp->derivatives().dy(2).transpose() * solution->v * cosa;
+                By_cos[i][2] += -msp->derivatives().dx(2).transpose() * solution->v * cosa;
+
+                Bx_sin[i][0] += msp->derivatives().dy(0).transpose() * solution->v * sina;
+                By_sin[i][0] += -msp->derivatives().dx(0).transpose() * solution->v * sina;
+                Bx_sin[i][1] += msp->derivatives().dy(1).transpose() * solution->v * sina;
+                By_sin[i][1] += -msp->derivatives().dx(1).transpose() * solution->v * sina;
+                Bx_sin[i][2] += msp->derivatives().dy(2).transpose() * solution->v * sina;
+                By_sin[i][2] += -msp->derivatives().dx(2).transpose() * solution->v * sina;
+            }
         }
 
         // Increment Position
@@ -187,5 +238,36 @@ void Model::simulate(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>  &fl
     }
     *position += Ns * dp;
 
-    conduction_losses = pp("ConductionLosses",Oe::VectorXd::Ones(previous[0].size()));
+#pragma omp parallel for
+    for (size_t i = 0; i < angle.size(); ++i) {
+        std::array<Eigen::VectorXd,3> B;
+        for (size_t j = 0; j < 3; ++j) {
+            Eigen::VectorXd B_cos = (Bx_cos[i][j].array() * Bx_cos[i][j].array() + By_cos[i][j].array() * By_cos[i][j].array()).sqrt();
+            Eigen::VectorXd B_sin = (Bx_sin[i][j].array() * Bx_sin[i][j].array() + By_sin[i][j].array() * By_sin[i][j].array()).sqrt();
+            B[j] = (B_cos.array() * B_cos.array() + B_sin.array() * B_sin.array()).sqrt();
+        }
+
+        for (auto cr : msp->domain()->regions()) {
+            auto r = cr.second;
+            if (r->material().magnetic_properties()->is_linear()) {
+                for (auto k : r->elements()) {
+                    for (size_t j = 0; j < 3; ++j) {
+                        B[j](k) = 0.0;
+                    }
+                }
+            }
+        }
+
+        //std::cout << B[0] << std::endl;
+
+        for (size_t j = 0; j < 3; ++j) {
+            losses[0](i) += (msp->weights()(j) * (99.89 * B[j].array().pow(1.809))).sum() * Params.ls * Params.Np * 1e-3;
+            losses[1](i) += (msp->weights()(j) * (3.788 * B[j].array().pow(1.500))).sum() * Params.ls * Params.Np * 1e-3;
+            losses[2](i) += (msp->weights()(j) * (0.301 * B[j].array().pow(2.000))).sum() * Params.ls * Params.Np * 1e-3;
+        }
+
+        losses[3](i) = pp("ConductionLosses", (current(i) * current(i) / Params.SlotFill / 2.0) * Oe::VectorXd::Ones(previous[0].size()));
+    }
+
+    ripple = (torque_max.array() - torque_min.array()) / 2.0;
 }

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, double_t &conduction_losses);
+    void simulate(Eigen::MatrixXd &torque, Eigen::MatrixXd &ripple, std::array<Eigen::MatrixXd,3>  &flux_cos, std::array<Eigen::MatrixXd,3>  &flux_sin, Eigen::MatrixXd angle, Eigen::MatrixXd current, std::array<Eigen::MatrixXd,4> &losses);
 
 protected:
     void create_sketch();

apps/Switched Winding Synchrel/parameters.cpp

 #include "parameters.h"
 
 Parameters::Parameters(Particle p) {
-    std::cerr << "Reminder: Implement analytical bridge thickness model" << std::endl;
+    // TODO: Fix meshing for small bridges (<0.5mm)
 
     Jmax = p["Jmax"].Position;
 
@@ -142,11 +142,5 @@ void Parameters::calculate_bridge_thickness() {
         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/torquespeedcurve.cpp

 #include "torquespeedcurve.h"
+#include <iostream>
 
-void torque_speed_curve(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>  &flux_cos, std::array<Eigen::MatrixXd,3>  &flux_sin, double_t DC_LINK, double_t Np, Eigen::ArrayXd &SPEED, Eigen::ArrayXd &TORQUE) {
+void torque_speed_curve(Eigen::MatrixXd &torque, Eigen::MatrixXd &ripple, std::array<Eigen::MatrixXd,4> &losses, std::array<Eigen::MatrixXd,3>  &flux_cos, std::array<Eigen::MatrixXd,3>  &flux_sin, double_t DC_LINK, double_t Np, Eigen::ArrayXd &SPEED, Eigen::ArrayXd &TORQUE, Eigen::ArrayXd &RIPPLE, Eigen::ArrayXd &LOSSES) {
     size_t Nr = (size_t)torque.rows();
     size_t Nc = (size_t)torque.cols();
 
@@ -18,7 +19,8 @@ void torque_speed_curve(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>
         }
     }
 
-    for (size_t s = 0; s != SPEED.size(); ++s) {
+#pragma omp parallel for
+    for (size_t s = 0; s < SPEED.size(); ++s) {
         Eigen::MatrixXd power = torque * (SPEED(s) * M_PI / 30.0);
 
         std::pair<size_t,size_t> ij;
@@ -49,7 +51,12 @@ void torque_speed_curve(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>
         Eigen::MatrixXd XY = Eigen::MatrixXd::Zero(9,9);
         Eigen::Vector3d xy; xy << -1.0, 0.0, 1.0;
         Eigen::VectorXd ztorque = Eigen::VectorXd::Zero(9);
+        Eigen::VectorXd zripple = Eigen::VectorXd::Zero(9);
         Eigen::VectorXd zspeed = Eigen::VectorXd::Zero(9);
+        Eigen::VectorXd zhys = Eigen::VectorXd::Zero(9);
+        Eigen::VectorXd zexc = Eigen::VectorXd::Zero(9);
+        Eigen::VectorXd zedd = Eigen::VectorXd::Zero(9);
+        Eigen::VectorXd zwnd = Eigen::VectorXd::Zero(9);
         for (size_t i = 0; i != 3; ++i) {
             size_t ii = i + ij.first - 1;
             for (size_t j = 0; j != 3; ++j) {
@@ -58,7 +65,12 @@ void torque_speed_curve(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>
                 size_t k = j + 3 * i;
 
                 ztorque(k) = torque(ii,jj);
+                zripple(k) = ripple(ii,jj);
                 zspeed(k) = max_speed(ii,jj);
+                zhys(k) = losses[0](ii,jj);
+                zexc(k) = losses[1](ii,jj);
+                zedd(k) = losses[2](ii,jj);
+                zwnd(k) = losses[3](ii,jj);
 
                 for (size_t m = 0; m != 3; ++m) {
                     for (size_t n = 0; n != 3; ++n) {
@@ -70,25 +82,39 @@ void torque_speed_curve(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>
         }
 
         Eigen::VectorXd torque_poly = XY.colPivHouseholderQr().solve(ztorque);
+        Eigen::VectorXd ripple_poly = XY.colPivHouseholderQr().solve(zripple);
         Eigen::VectorXd speed_poly = XY.colPivHouseholderQr().solve(zspeed);
+        Eigen::VectorXd chys_poly = XY.colPivHouseholderQr().solve(zhys);
+        Eigen::VectorXd cexc_poly = XY.colPivHouseholderQr().solve(zexc);
+        Eigen::VectorXd cedd_poly = XY.colPivHouseholderQr().solve(zedd);
+        Eigen::VectorXd cwnd_poly = XY.colPivHouseholderQr().solve(zwnd);
 
         for (size_t i = 0; i != 201; ++i) {
             double_t x = (i-100.0) / 100.0;
             for (size_t j = 0; j != 201; ++j) {
                 double_t y = (j-100.0) / 100.0;
 
-                double_t torque{0.0}, speed{0.0};
+                double_t torque{0.0}, ripple{0.0}, speed{0.0}, chys{0.0}, cexc{0.0}, cedd{0.0}, cwnd{0.0};
                 for (size_t m = 0; m != 3; ++m) {
                     for (size_t n = 0; n != 3; ++n) {
                         size_t p = n + 3 * m;
 
                         torque += torque_poly(p) * pow(x,m) * pow(y,n);
+                        ripple += ripple_poly(p) * pow(x,m) * pow(y,n);
                         speed += speed_poly(p) * pow(x,m) * pow(y,n);
+                        chys += chys_poly(p) * pow(x,m) * pow(y,n);
+                        cexc += cexc_poly(p) * pow(x,m) * pow(y,n);
+                        cedd += cedd_poly(p) * pow(x,m) * pow(y,n);
+                        cwnd += cwnd_poly(p) * pow(x,m) * pow(y,n);
                     }
                 }
 
                 if (speed >= SPEED[s] && torque >= TORQUE[s]) {
                     TORQUE[s] = torque;
+                    RIPPLE[s] = ripple;
+
+                    double_t fe = SPEED[s] * 1e3 / 60.0 * Np / 2.0;
+                    LOSSES[s] = chys * fe + cexc * std::pow(fe,1.5) + cedd * std::pow(fe,2.0) + cwnd;
                 }
             }
         }

apps/Switched Winding Synchrel/torquespeedcurve.h

@@ -3,6 +3,6 @@
 
 #include "Eigen"
 
-void torque_speed_curve(Eigen::MatrixXd &torque, std::array<Eigen::MatrixXd,3>  &flux_cos, std::array<Eigen::MatrixXd,3>  &flux_sin, double_t DC_LINK, double_t Np, Eigen::ArrayXd &SPEED, Eigen::ArrayXd &TORQUE);
+void torque_speed_curve(Eigen::MatrixXd &torque, Eigen::MatrixXd &ripple, std::array<Eigen::MatrixXd,4> &losses, std::array<Eigen::MatrixXd,3>  &flux_cos, std::array<Eigen::MatrixXd,3>  &flux_sin, double_t DC_LINK, double_t Np, Eigen::ArrayXd &SPEED, Eigen::ArrayXd &TORQUE, Eigen::ArrayXd &RIPPLE, Eigen::ArrayXd &LOSSES);
 
 #endif //OERSTED_TORQUESPEEDCURVE_H

build.sh

@@ -3,16 +3,16 @@ cd build
 
 cmake .. && make
 
-./test/run_tests
+#./test/run_tests
 
-mkdir -p coverage
+#mkdir -p coverage
 
-COV_FILE="coverage/coverage.info"
-COV_DIR=".coverage"
-COV_SITE=$COV_DIR/"index.html"
-RM_PREFIX="${PWD%/*}/src"
+#COV_FILE="coverage/coverage.info"
+#COV_DIR=".coverage"
+#COV_SITE=$COV_DIR/"index.html"
+#RM_PREFIX="${PWD%/*}/src"
 
-lcov -c -d .. -o $COV_FILE
-lcov -r $COV_FILE "*Oersted/lib/*" "*Oersted/test/*" "/usr/include/*" "/usr/lib/*" -o $COV_FILE
-genhtml $COV_FILE -o $COV_DIR -p $RM_PREFIX
-google-chrome $COV_SITE
+#lcov -c -d .. -o $COV_FILE
+#lcov -r $COV_FILE "*Oersted/lib/*" "*Oersted/test/*" "/usr/include/*" "/usr/lib/*" -o $COV_FILE
+#genhtml $COV_FILE -o $COV_DIR -p $RM_PREFIX
+#google-chrome $COV_SITE

src/Elements/CMakeLists.txt

@@ -9,6 +9,6 @@ set(SOURCE_FILES
         ./src/Element.cpp   ./src/Element.h
         ./src/Facet.cpp     ./src/Facet.h)
 
-add_library(elements SHARED ${SOURCE_FILES})
+add_library(elements STATIC ${SOURCE_FILES})
 
 target_link_libraries(elements)
\ No newline at end of file

src/Geometry/CMakeLists.txt

@@ -6,6 +6,6 @@ set(SOURCE_FILES
         ./src/Ray.cpp ./src/Ray.h
         ./src/Circle.cpp ./src/Circle.h)
 
-add_library(geometry SHARED ${SOURCE_FILES})
+add_library(geometry STATIC ${SOURCE_FILES})
 
 target_link_libraries(geometry)
\ No newline at end of file

src/Materials/CMakeLists.txt

@@ -8,6 +8,6 @@ set(SOURCE_FILES
         ./src/MaterialProperties.cpp ./src/MaterialProperties.h
         )
 
-add_library(materials SHARED ${SOURCE_FILES})
+add_library(materials STATIC ${SOURCE_FILES})
 
 target_link_libraries(materials)
\ No newline at end of file

src/Materials/src/MaterialProperties.h

@@ -11,7 +11,12 @@
 
 #include "PhysicsConstants.h"
 
-class MagneticMaterialInterface {
+class MaterialInterface {
+public:
+    virtual bool is_linear() const = 0;
+};
+
+class MagneticMaterialInterface : public MaterialInterface {
 public:
     virtual ~MagneticMaterialInterface(){};
 
@@ -24,6 +29,8 @@ class LinearIsotropicMagneticMaterial : public MagneticMaterialInterface {
 public:
     LinearIsotropicMagneticMaterial(double_t relative_permeability) : Reluctivity(1.0 / (mu_0 * relative_permeability)) {};
 
+    bool is_linear() const override { return true; };
+
     void h_from_b(std::vector<size_t> const &index, Oe::ArrayXd &Fx, Oe::ArrayXd &Fy, Oe::ArrayXd &dFxdx, Oe::ArrayXd &dFydy, Oe::ArrayXd &dFxdy) override {
         for(size_t const &i : index) {
             Fx(i) *= Reluctivity;
@@ -42,7 +49,7 @@ public:
     }
 
 protected:
-    double Reluctivity;
+    double_t Reluctivity;
 };
 
 class MaterialProperties {
@@ -50,6 +57,8 @@ public:
     // TODO: Add std::string Name; property + a constructor which loads data saved in a file on disk
     MaterialProperties(std::shared_ptr<MagneticMaterialInterface> magnetic) : MagneticProperties{magnetic} {};
 
+    std::shared_ptr<MagneticMaterialInterface> magnetic_properties() const { return MagneticProperties; };
+
     template<typename... Args>
     void h_from_b(Args&&... args) { return MagneticProperties->h_from_b(std::forward<Args>(args)...); }
 
@@ -83,6 +92,8 @@ public:
         MSaturation = (s[0] * BSaturation + s[1]) * BSaturation + s[2];
     }
 
+    bool is_linear() const override { return false; };
+
     void h_from_b(std::vector<size_t> const &index, Oe::ArrayXd &Fx, Oe::ArrayXd &Fy, Oe::ArrayXd &dFxdx, Oe::ArrayXd &dFydy, Oe::ArrayXd &dFxdy) override {
         for(size_t const &i : index) {
             const double_t &Bx = Fx(i);

src/Matrix/CMakeLists.txt

@@ -13,6 +13,6 @@ set(SOURCE_FILES
         ./src/SecondDerivativeMatrixGroup.cpp ./src/SecondDerivativeMatrixGroup.h
         )
 
-add_library(matrix SHARED ${SOURCE_FILES})
+add_library(matrix STATIC ${SOURCE_FILES})
 
 target_link_libraries(matrix)
\ No newline at end of file

src/Mesh/CMakeLists.txt

@@ -11,6 +11,6 @@ set(SOURCE_FILES
         ./src/Mesh.h                ./src/Mesh.cpp
         ./src/DartTriangle.h        ./src/DartTriangle.cpp)
 
-add_library(mesh SHARED ${SOURCE_FILES})
+add_library(mesh STATIC ${SOURCE_FILES})
 
-target_link_libraries(mesh)
\ No newline at end of file
+target_link_libraries(mesh geometry elements)
\ No newline at end of file

src/Model Templates/CMakeLists.txt

@@ -9,6 +9,6 @@ set(SOURCE_FILES
         ./src/PolarModelTemplate.cpp ./src/PolarModelTemplate.h
         )
 
-add_library(model_templates SHARED ${SOURCE_FILES})
+add_library(model_templates STATIC ${SOURCE_FILES})
 
-target_link_libraries(model_templates)
\ No newline at end of file
+target_link_libraries(model_templates sketch)
\ No newline at end of file

src/Optimization/CMakeLists.txt

@@ -10,6 +10,6 @@ set(SOURCE_FILES