Further refactoring of Physics to move some dependencies to the FiniteElementMesh class

Changes the way number of elements and unknowns are calculated in Physics
    This still needs some work
Adds FunctionArgument class for passing arguments to Forcing objects with multiple parameters
Adds test demonstrating parallelism using OpenMP

CMakeLists.txt

@@ -3,8 +3,9 @@ cmake_minimum_required(VERSION 3.2)
 project(Oersted)
 
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=gnu++17")
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O0 --coverage")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O3 --coverage")
 #set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp")
 
 include_directories(./lib/)
 include_directories(./lib/Eigen/)

src/Physics/src/FiniteElementMesh.h

@@ -23,8 +23,6 @@ public:
 
     virtual ~FiniteElementMeshInterface() {};
 
-    size_t size_nodes() const { return Nodes.size(); };
-
     auto &boundaries() { return Boundaries; };
 
     auto &boundary(size_t i) { return Boundaries[i]; };
@@ -75,13 +73,11 @@ public:
         }
     }
 
-    size_t nodes_size() const { return Nodes.size(); };
+    DiagonalMatrixGroup diagonal_matrix_group() const { return DiagonalMatrixGroup{quadrature_size(), elements_size()}; };
 
-    virtual size_t elements_size() const = 0;
+    SparseMatrixGroup sparse_matrix_group() const { return SparseMatrixGroup{quadrature_size(), nodes_size(), elements_size(), nodes_per_element()}; };
 
-    virtual void write_scalar(Eigen::VectorXd const &v, std::string path, std::string file_name) const = 0;
-
-    virtual void write_vector(Eigen::ArrayXd const &Fx, Eigen::ArrayXd const &Fy, std::string path, std::string file_name) const = 0;
+    DerivativeMatrixGroup derivative_matrix_group() const { return DerivativeMatrixGroup{quadrature_size(), nodes_size(), elements_size(), nodes_per_element()}; };
 
     virtual DiagonalMatrixGroup determinant() const = 0;
 
@@ -91,6 +87,10 @@ public:
 
     virtual std::vector<std::vector<XY>> quadrature_points() const = 0;
 
+    size_t nodes_size() const { return Nodes.size(); };
+
+    virtual size_t elements_size() const = 0;
+
     virtual size_t element_order() const = 0;
 
     virtual size_t nodes_per_element() const = 0;
@@ -101,6 +101,10 @@ public:
 
     virtual double_t quadrature_weight(size_t i) const = 0;
 
+    virtual void write_scalar(Eigen::VectorXd const &v, std::string path, std::string file_name) const = 0;
+
+    virtual void write_vector(Eigen::ArrayXd const &Fx, Eigen::ArrayXd const &Fy, std::string path, std::string file_name) const = 0;
+
 protected:
     std::vector<XY> Nodes;
 

src/Physics/src/Forcing.h

@@ -3,6 +3,7 @@
 
 #include <functional>
 #include <vector>
+#include <unordered_map>
 
 #include "Eigen"
 
@@ -12,16 +13,26 @@
 #include "MatrixGroup.h"
 #include "FiniteElementMesh.h"
 
+class FunctionArguments {
+public:
+    void add_argument(std::string key, double_t value) { Arguments.emplace(key, value); };
+
+    double_t &operator[](std::string key) { return Arguments[key]; };
+
+protected:
+    std::unordered_map<std::string, double_t> Arguments;
+};
+
 class Forcing {
 public:
     virtual ~Forcing() {};
 
-    virtual Eigen::VectorXd operator()(double const t) const = 0;
+    virtual Eigen::VectorXd operator()(FunctionArguments fargs) const = 0;
 };
 
 class HomogeneousForcing : public Forcing {
 public:
-    HomogeneousForcing(std::function<double(double)> f,
+    HomogeneousForcing(std::function<double(FunctionArguments)> f,
                        std::vector<std::shared_ptr<DiscreteRegion>> regions,
                        SparseMatrixGroup const &basis,
                        DiagonalMatrixGroup const &weights)
@@ -37,20 +48,20 @@ public:
         }
     };
 
-    Eigen::VectorXd operator()(double const t) const override {
+    Eigen::VectorXd operator()(FunctionArguments fargs) const override {
         Eigen::VectorXd output = (Basis[0] * (Weights(0).matrix().asDiagonal() * Indicator)); // TODO: ?SparseVector?
         for (size_t i = 1; i != Basis.size(); ++i) {
             output += (Basis[i] * (Weights(i).matrix().asDiagonal() * Indicator));
         }
 
-        output *= Function(t);
+        output *= Function(fargs);
         return output;
     };
 
     auto const &indicator() const { return Indicator; };
 
 protected:
-    std::function<double(double)> Function;
+    std::function<double(FunctionArguments)> Function;
 
     std::vector<std::shared_ptr<DiscreteRegion>> Regions;
 

src/Physics/src/MatrixGroup.h

@@ -8,6 +8,8 @@
 
 class SparseMatrixGroup {
 public:
+    SparseMatrixGroup() {};
+
     SparseMatrixGroup(size_t const Q, size_t const rows, size_t const cols, size_t const nnz) {
         Matrices.resize(Q);
         for (size_t i = 0; i != Matrices.size(); ++i) {
@@ -36,6 +38,8 @@ protected:
 
 class DiagonalMatrixGroup {
 public:
+    DiagonalMatrixGroup() {};
+
     DiagonalMatrixGroup(size_t Q, size_t const dim) {
         Matrices.resize(Q);
         for (size_t i = 0; i != Matrices.size(); ++i) {
@@ -55,6 +59,8 @@ protected:
 
 class DerivativeMatrixGroup {
 public:
+    DerivativeMatrixGroup() {};
+
     DerivativeMatrixGroup(size_t const Q, size_t const rows, size_t const cols, size_t const nnz) : Dx{Q, rows, cols, nnz}, Dy{Q, rows, cols, nnz} {};
 
     auto const &dx(size_t const q) const { return Dx[q]; };

src/Physics/src/Physics.h

@@ -29,21 +29,12 @@ public:
     virtual void apply_conditions() = 0;
      */
 
-    PhysicsInterface(std::shared_ptr<FiniteElementMeshInterface> d) : Domain{d},
-                                                                      ElementWeights{d->quadrature_size(), d->elements_size()},
-                                                                      Derivatives{d->quadrature_size(), d->nodes_size(), d->elements_size(), d->nodes_per_element()},
-                                                                      Basis{d->quadrature_size(), d->nodes_size(), d->elements_size(), d->nodes_per_element()} {
+    PhysicsInterface(std::shared_ptr<FiniteElementMeshInterface> d) : Domain{d} {
         std::cout << "TODO: Think about storing matrices in the FiniteElementMesh object since this class is simply forwarding it's properties" << std::endl;
         std::cout << "TODO: Rather, give FiniteElementMesh methods which produce ElementWeight, Derivative, and Basis matrices of appropriate sizes" << std::endl;
         std::cout << "TODO: That is, FiniteElementMesh is essentially a MatrixFactory" << std::endl;
     };
 
-    double time() const { return Time; };
-
-    double &time() { return Time; };
-
-    void time(double t) { Time = t; };
-
     std::shared_ptr<FiniteElementMeshInterface> domain() const { return Domain; };
 
     DiagonalMatrixGroup const &weights() { return ElementWeights; };
@@ -54,8 +45,12 @@ public:
 
     Eigen::VectorXd recover_boundary(Eigen::VectorXd const &v) const { return BoundaryConditionMatrix.transpose() * v; };
 
+    virtual size_t unknowns_size() const = 0;
+
+    virtual size_t elements_size() const = 0;
+
 protected:
-    double Time{0.0};
+    size_t Unknowns{0};
 
     std::shared_ptr<FiniteElementMeshInterface> Domain;
 
@@ -71,13 +66,9 @@ protected:
 
 class Magnetostatic : public PhysicsInterface {
 public:
-    Magnetostatic(std::shared_ptr<FiniteElementMeshInterface> d) : PhysicsInterface{d}, Unknowns{d->nodes_size()}, Elements{d->elements_size()} {}
-
-    Eigen::SparseMatrix<double> init_unknown_matrix() const { return Eigen::SparseMatrix<double>(Unknowns, Unknowns); }
+    Magnetostatic(std::shared_ptr<FiniteElementMeshInterface> d) : PhysicsInterface{d} { build_matrices(); };
 
-    Eigen::DiagonalMatrix<double, Eigen::Dynamic> init_element_matrix() const { return Eigen::DiagonalMatrix<double, Eigen::Dynamic>(Elements); }
-
-    std::shared_ptr<Solution> initialize() {
+    std::shared_ptr<Solution> new_solution() {
         auto s = std::make_shared<Solution>();
 
         s->J = init_unknown_matrix(); // Jacobian
@@ -103,20 +94,20 @@ public:
         apply_conditions();
     }
 
-    Eigen::VectorXd init_unknown_vector() const { return Eigen::VectorXd(Unknowns); }
+    Eigen::SparseMatrix<double> init_unknown_matrix() const { return Eigen::SparseMatrix<double>(unknowns_size(), unknowns_size()); }
 
-    Eigen::VectorXd init_element_vector() const { return Eigen::VectorXd(Elements); }
+    Eigen::VectorXd init_unknown_vector() const { return Eigen::VectorXd(unknowns_size()); }
 
-    Eigen::ArrayXd init_unknown_array() const { return Eigen::ArrayXd(Unknowns); }
+    Eigen::VectorXd init_element_vector() const { return Eigen::VectorXd(elements_size()); }
 
-    Eigen::ArrayXd init_element_array() const { return Eigen::ArrayXd(Elements); }
+    Eigen::ArrayXd init_element_array() const { return Eigen::ArrayXd(elements_size()); }
 
-    void solve(std::shared_ptr<Solution> s) {
+    void solve(std::shared_ptr<Solution> s, FunctionArguments fargs) {
         // TODO: Template on solver type
         // TODO: Save J,v,r,f,Fx,Fy,dFxdx,dFydy,dFxdy to reduce number of memory allocations?
         // TODO: This method is fairly general
 
-        calculate_forcing(s->f);
+        calculate_forcing(s->f, fargs);
 
 
         linearize(s->J, s->v, s->r, s->f, s->Fx, s->Fy, s->dFxdGx, s->dFydGy, s->dFxdGy);
@@ -150,7 +141,7 @@ public:
         r = -f;
         J.setZero();
         for (size_t i = 0; i != Domain->quadrature_size(); ++i) {
-            // Caclulate flux density
+            // Calculate flux density
             Fx = this->Derivatives.dy(i).transpose() * v;
             Fy = -this->Derivatives.dx(i).transpose() * v;
 
@@ -163,12 +154,11 @@ public:
             r += this->Derivatives.dy(i) * (this->ElementWeights(i) * Fx).matrix()
                  - this->Derivatives.dx(i) * (this->ElementWeights(i) * Fy).matrix(); // note: weak form introduces a negative sign into double-curl operator
 
-            // Calculate jacobian
+            // Calculate Jacobian
             J += this->Derivatives.dx(i) * (this->ElementWeights(i) * dFyGy).matrix().asDiagonal() * this->Derivatives.dx(i).transpose()
                  + this->Derivatives.dy(i) * (this->ElementWeights(i) * dFxGx).matrix().asDiagonal() * this->Derivatives.dy(i).transpose()
                  - this->Derivatives.dx(i) * (this->ElementWeights(i) * dFxGy).matrix().asDiagonal() * this->Derivatives.dy(i).transpose()
                  - this->Derivatives.dy(i) * (this->ElementWeights(i) * dFxGy).matrix().asDiagonal() * this->Derivatives.dx(i).transpose();
-
         }
     }
 
@@ -180,16 +170,19 @@ public:
         for (size_t i = 0; i != Domain->quadrature_size(); ++i) {
             this->ElementWeights(i) *= Domain->quadrature_weight(i);//TriangleQuadrature<(Domain->quadrature_order())>::w[i];
         }
+
+        this->BoundaryConditionMatrix.resize(this->Domain->nodes_size(), this->Domain->nodes_size());
+        this->BoundaryConditionMatrix.setIdentity();
     }
 
-    void calculate_forcing(Eigen::VectorXd &f) {
+    void calculate_forcing(Eigen::VectorXd &f, FunctionArguments fargs) {
         f.setZero();
         for (size_t i = 0; i != this->ForcingCondtions.size(); ++i) {
-            f += this->ForcingCondtions[i]->operator()(this->Time);
+            f += this->ForcingCondtions[i]->operator()(fargs);
         }
     }
 
-    void add_current_density(std::function<double(double)> func, std::vector<std::shared_ptr<Contour const>> contours) {
+    void add_current_density(std::function<double(FunctionArguments)> func, std::vector<std::shared_ptr<Contour const>> contours) {
         std::vector<std::shared_ptr<DiscreteRegion>> regions;
         for (auto &c : contours) {
             regions.push_back(this->Domain->region(c));
@@ -198,7 +191,7 @@ public:
         this->ForcingCondtions.push_back(std::make_shared<HomogeneousForcing>(func, regions, this->Basis, this->ElementWeights));
     }
 
-    void add_current_density(std::function<double(double)> func, std::vector<std::shared_ptr<DiscreteRegion>> regions) {
+    void add_current_density(std::function<double(FunctionArguments)> func, std::vector<std::shared_ptr<DiscreteRegion>> regions) {
         this->ForcingCondtions.push_back(std::make_shared<HomogeneousForcing>(func, regions, this->Basis, this->ElementWeights));
     }
 
@@ -227,8 +220,7 @@ public:
         return pbc;
     }
 
-    auto add_periodic_boundary(std::vector<std::shared_ptr<DiscreteBoundary>> b0, std::vector<std::shared_ptr<DiscreteBoundary>> b1, std::vector<bool> orientation,
-                               bool antiperiodic) {
+    auto add_periodic_boundary(std::vector<std::shared_ptr<DiscreteBoundary>> b0, std::vector<std::shared_ptr<DiscreteBoundary>> b1, std::vector<bool> orientation, bool antiperiodic) {
         auto pbc = std::make_shared<PeriodicBoundaryCondition>(b0, b1, orientation, antiperiodic);
         this->BoundaryConditions.push_back(pbc);
         return pbc;
@@ -252,9 +244,6 @@ public:
         // TODO: This is a general method that should be in a superclass
 
         // Apply permutation matrices
-        this->BoundaryConditionMatrix.resize(this->Domain->nodes_size(), this->Domain->nodes_size());
-        this->BoundaryConditionMatrix.setIdentity();
-
         for (size_t i = 0; i != this->BoundaryConditions.size(); ++i) {
             this->BoundaryConditions[i]->apply(this->BoundaryConditionMatrix);
         }
@@ -282,21 +271,15 @@ public:
 
         // Reduce boundary conditions
         this->BoundaryConditionMatrix = reduction_matrix * this->BoundaryConditionMatrix;
-
         this->Basis.transform(this->BoundaryConditionMatrix);
         this->Derivatives.transform(this->BoundaryConditionMatrix);
-        this->Unknowns = this->Domain->nodes_size() - reduction_index.size();
     }
 
-    size_t unknowns() const { return Unknowns; }
+    size_t unknowns_size() const override { return (size_t)BoundaryConditionMatrix.rows(); }
 
-    size_t elements() const { return Elements; }
+    size_t elements_size() const override { return Domain->elements_size(); }
 
 protected:
-
-private:
-    size_t Unknowns;
-    size_t Elements;
 };
 
 #endif //OERSTED_PHYSICS_H

test/LibraryIntegration/Library_Integration_Test.hpp

@@ -8,4 +8,6 @@
 #include "Quadrature.hpp"
 #include "Physics.hpp"
 
+#include <omp.h>
+
 #endif //OERSTED_TEST_LIBRARY_INTEGRATION_H

test/LibraryIntegration/Mesh_To_FEM_Test.cpp

@@ -69,7 +69,9 @@ TEST(Full_Circle, Magnetostatic_Uniform_Current_Density) {
     // Create magnetostatic physics
     Magnetostatic msph{fem};
 
-    msph.add_current_density([](double t) { return 1.0 / (2.0 * M_PI * 1e-7); }, {fem->region(0)});
+    FunctionArguments fargs;
+
+    msph.add_current_density([](FunctionArguments args) { 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)});
 
@@ -90,8 +92,7 @@ TEST(Full_Circle, Magnetostatic_Uniform_Current_Density) {
 
     //
     // Set time
-    msph.time(0.0);
-    msph.calculate_forcing(f);
+    msph.calculate_forcing(f, fargs);
 
     // Linearize
     v.setZero();
@@ -332,7 +333,9 @@ TEST_F(Sixth_Circle_Mirror_Copy, Magnetostatic_Uniform_Current_Density_Periodic)
     // Create Physics
     Magnetostatic msph{fem};
 
-    msph.add_current_density([](double t) { return 1.0 / (2.0 * M_PI * 1e-7); }, {current_density_contour});
+    FunctionArguments fargs;
+
+    msph.add_current_density([](FunctionArguments args) { return 1.0 / (2.0 * M_PI * 1e-7); }, {current_density_contour});
 
     msph.add_magnetic_insulation(outer_boundary);
 
@@ -354,7 +357,7 @@ TEST_F(Sixth_Circle_Mirror_Copy, Magnetostatic_Uniform_Current_Density_Periodic)
 
     v.setZero();
 
-    msph.calculate_forcing(f);
+    msph.calculate_forcing(f, fargs);
 
     msph.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
 
@@ -430,7 +433,9 @@ TEST_F(Sixth_Circle_Mirror_Copy, Magnetostatic_Uniform_Current_Density_Antieriod
     // Create Physics
     Magnetostatic msph{fem};
 
-    msph.add_current_density([](double t) { return 1.0 / (2.0 * M_PI * 1e-7); }, {current_density_contour});
+    FunctionArguments fargs;
+
+    msph.add_current_density([](FunctionArguments args) { return 1.0 / (2.0 * M_PI * 1e-7); }, {current_density_contour});
 
     msph.add_magnetic_insulation(outer_boundary);
 
@@ -452,7 +457,7 @@ TEST_F(Sixth_Circle_Mirror_Copy, Magnetostatic_Uniform_Current_Density_Antieriod
 
     v.setZero();
 
-    msph.calculate_forcing(f);
+    msph.calculate_forcing(f, fargs);
 
     msph.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
 
@@ -706,7 +711,9 @@ TEST_F(Sixth_Circle, Magnetostatic_Uniform_Current_Density_Periodic_Rotation) {
         }
     }
 
-    msph.add_current_density([](double t) { return 1.0 / (2.0 * M_PI * 1e-7); }, {current_density_contour});
+    FunctionArguments fargs;
+
+    msph.add_current_density([](FunctionArguments args) { return 1.0 / (2.0 * M_PI * 1e-7); }, {current_density_contour});
 
     msph.add_magnetic_insulation(outer_boundary);
 
@@ -729,7 +736,7 @@ TEST_F(Sixth_Circle, Magnetostatic_Uniform_Current_Density_Periodic_Rotation) {
 
         v.setZero();
 
-        msph.calculate_forcing(f);
+        msph.calculate_forcing(f, fargs);
 
         msph.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
 
@@ -812,7 +819,9 @@ TEST_F(Sixth_Circle, Magnetostatic_Uniform_Current_Density_Antiperiodic_Rotation
     double_t angle = 0.0;
     double_t delta_angle = (2 * M_PI) / (6.0 * position->size());
 
-    msph.add_current_density([](double t) { return 1.0 / (2.0 * M_PI * 1e-7); }, {current_density_contour});
+    FunctionArguments fargs;
+
+    msph.add_current_density([](FunctionArguments args) { return 1.0 / (2.0 * M_PI * 1e-7); }, {current_density_contour});
 
     msph.add_magnetic_insulation(outer_boundary);
 
@@ -835,7 +844,7 @@ TEST_F(Sixth_Circle, Magnetostatic_Uniform_Current_Density_Antiperiodic_Rotation
 
         v.setZero();
 
-        msph.calculate_forcing(f);
+        msph.calculate_forcing(f, fargs);
 
         msph.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
 
@@ -865,7 +874,7 @@ TEST_F(Sixth_Circle, Magnetostatic_Uniform_Current_Density_Antiperiodic_Rotation
 
         double_t r_max = -(2.0 * b_f) / (3.0 * a_f);
         double_t tol_a = (a_f * r_max + b_f) * r_max * r_max * 0.02;
-        for (size_t i = 0; i != fem->size_nodes(); ++i) {
+        for (size_t i = 0; i != fem->nodes_size(); ++i) {
             XY const &n = fem->node(i);
             double_t r = std::hypot(n.x(), n.y());
             double_t a = std::atan2(n.y(), n.x());
@@ -1148,31 +1157,32 @@ TEST_F(Salient_Pole_Synchrel, Test) {
     Magnetostatic msph{fem};
 
     // Set material properties
-    //MaterialProperties electrical_steel = MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1000.0));
-    //MaterialProperties electrical_steel = MaterialProperties(std::make_shared<NonlinearIsotropicMagneticMaterial>());
     MaterialProperties electrical_steel = JFE_35JN210();
-
     fem->region(rotor_iron)->material(electrical_steel);
     fem->region(stator_iron)->material(electrical_steel);
 
     // Conditions
-    msph.add_current_density([](double t) { return -2e6 * sin(2 * M_PI * t - M_PI * 2.0 / 3.0); }, {phase_a_contour});
-    msph.add_current_density([](double t) { return 2e6 * sin(2 * M_PI * t); }, {phase_b_contour});
-    msph.add_current_density([](double t) { return -2e6 * sin(2 * M_PI * t + M_PI * 2.0 / 3.0); }, {phase_c_contour});
+    FunctionArguments fargs;
+    fargs.add_argument("t",0.0);
+    fargs.add_argument("J",2e6);
+
+    msph.add_current_density([](FunctionArguments args) { return -args["J"] * sin(2 * M_PI * args["t"] - M_PI * 2.0 / 3.0); }, {phase_a_contour});
+    msph.add_current_density([](FunctionArguments args) { return +args["J"] * sin(2 * M_PI * args["t"] + M_PI * 0.0 / 3.0); }, {phase_b_contour});
+    msph.add_current_density([](FunctionArguments args) { return -args["J"] * sin(2 * M_PI * args["t"] + M_PI * 2.0 / 3.0); }, {phase_c_contour});
 
     msph.add_magnetic_insulation(outer_boundary);
 
     msph.add_periodic_boundary(periodic_boundary, true);
 
     auto position = msph.add_sliding_interface(continuous_airgap, true);
+    double_t t = 0.0;
     double_t dt = 1.0 / (2.0 * position->size());
 
-    msph.assemble();
-
     // Solve
-    auto solution = msph.initialize();
+    msph.assemble();
     for (size_t iter = 0; iter != position->size() / 32; ++iter) {
-        msph.solve(solution);
+        auto solution = msph.new_solution(); // TODO: What about warm start using old initial condition?
+        msph.solve(solution, fargs);
 
         //Verify equation is solved
         EXPECT_LE(solution->r.norm(), FLT_EPSILON * solution->f.norm() * sqrt(solution->f.size()));
@@ -1187,8 +1197,71 @@ TEST_F(Salient_Pole_Synchrel, Test) {
 
         // Increment Position
         *position += 32;
-        msph.time() += 32 * dt;
+        t += 32 * dt;
 
         msph.assemble();
     }
+}
+
+TEST_F(Salient_Pole_Synchrel, Test_OpenMP) {
+    #ifndef _OPENMP
+    std::cout << "Test must be compiled with -fopenmp" << std::endl;
+    FAIL();
+    #endif
+
+    omp_set_num_threads(omp_get_max_threads() / 4);
+    size_t N = (size_t)(std::sqrt(omp_get_max_threads() / 4));
+
+    // Input Arguments
+    std::vector <double_t> t;
+    std::vector <double_t> J;
+    for (size_t i = 0; i != N; ++i) {
+        for (size_t j = 0; j != N; ++j) {
+            t.push_back(i / (N - 1.0));
+            J.push_back(2.0e6 * j / (N - 1.0));
+        }
+    }
+
+    Magnetostatic msph{fem};
+
+    // Materials
+    MaterialProperties electrical_steel = JFE_35JN210();
+    fem->region(rotor_iron)->material(electrical_steel);
+    fem->region(stator_iron)->material(electrical_steel);
+
+    // Conditions
+    msph.add_current_density([](FunctionArguments args) { return -args["J"] * sin(2 * M_PI * args["t"] - M_PI * 2.0 / 3.0); }, {phase_a_contour});
+    msph.add_current_density([](FunctionArguments args) { return +args["J"] * sin(2 * M_PI * args["t"] + M_PI * 0.0 / 3.0); }, {phase_b_contour});
+    msph.add_current_density([](FunctionArguments args) { return -args["J"] * sin(2 * M_PI * args["t"] + M_PI * 2.0 / 3.0); }, {phase_c_contour});
+
+    msph.add_magnetic_insulation(outer_boundary);
+
+    msph.add_periodic_boundary(periodic_boundary, true);
+
+    auto position = msph.add_sliding_interface(continuous_airgap, true);
+    double_t dt = 1.0 / (2.0 * position->size());
+
+    // Solve
+    for (size_t iter = 0; iter != position->size() / 32; ++iter) {
+        std::cout << iter << std::endl;
+        msph.assemble(); // assemble should not be called within an OpenMP parallel construct since Physics objects are global
+
+        #pragma omp parallel for
+        for (size_t i = 0; i < t.size(); ++i) {
+            FunctionArguments fargs;
+            fargs.add_argument("t", t[i]);
+            fargs.add_argument("J", J[i]);
+
+            auto solution = msph.new_solution();
+            // TODO: What about warm start using previous solution as initial condition? #pragma omp parallel for schedule(static,1) maybe not
+            // TODO: Need to come up with an interleaving strategy based on thread id
+            msph.solve(solution, fargs);
+
+            //Verify equation is solved
+            EXPECT_LE(solution->r.norm(), FLT_EPSILON * solution->f.norm() * sqrt(solution->f.size()));
+        }
+
+        // Increment Position
+        *position += 32;
+    }
 }
\ No newline at end of file

test/Physics/Finite_Element_Mesh_Test.cpp

@@ -75,9 +75,11 @@ TEST_F(MasterTriangle, magnetostatic) {
     auto femesh = get_mesh<2>();
 
     Magnetostatic ms = Magnetostatic(femesh);
-    ms.time(1.0);
 
-    auto ff = [](double t) { return 1.0 * t; };
+    FunctionArguments fargs;
+    fargs.add_argument("t", 1.0);
+
+    auto ff = [](FunctionArguments args) { return 1.0 * args["t"]; };
 
     ms.add_current_density(ff, reg);
 
@@ -93,10 +95,10 @@ TEST_F(MasterTriangle, magnetostatic) {
     auto dFydy = ms.init_element_array();
     auto dFxdy = ms.init_element_array();
 
-    ms.calculate_forcing(f);
+    ms.calculate_forcing(f, fargs);
     ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
 
-    PostProcessor mspp{std::make_shared<PhysicsInterface>(ms)};
+    PostProcessor mspp{std::make_shared<Magnetostatic>(ms)};
 
     //auto ff = [](double t) { return 1.0 * t; };
 
@@ -370,14 +372,15 @@ TEST_F(SimpleSquare, magnetostatic_forcing) {
     auto femesh = get_mesh<2>();
 
     Magnetostatic ms{femesh};
-    ms.time(1.0);
     ms.build_matrices();
 
-    auto ff = [](double t) { return 1.0 * t; };
+    FunctionArguments fargs;
+    fargs.add_argument("t",1.0);
+    auto ff = [](FunctionArguments args) { return 1.0 * args["t"]; };
     auto f = ms.init_unknown_vector();
 
     ms.add_current_density(ff, {reg[0]});
-    ms.calculate_forcing(f);
+    ms.calculate_forcing(f, fargs);
     EXPECT_NEAR(f(0), 1.0 / 6.0, FLT_EPSILON);
     EXPECT_NEAR(f(1), 1.0 / 6.0, FLT_EPSILON);
     EXPECT_NEAR(f(2), 1.0 / 6.0, FLT_EPSILON);