Implements SpatialForcing which allows, e.g., spatially varying current...

Implements SpatialForcing which allows, e.g., spatially varying current densities defined by lambda functions (e.g. sinusoids)
Implements tests for SpatialForcing using sinusoidal functions in test_Square and test_Circle
Begins implementation of tests on two circular annuli (test_Annuli)

src/Matrix/src/DiagonalMatrixGroup.h

@@ -18,7 +18,7 @@ public:
 
     auto &operator()(size_t const q) { return Matrices[q]; };
 
-    size_t size() const { return (size_t)(Matrices[0].size()); };
+    size_t size() const { return (size_t)(Matrices[0].size()); }; // TODO: This should return Matrices.size() for consistency
 
 protected:
     std::vector<Oe::ArrayXd> Matrices;

src/Physics/CMakeLists.txt

@@ -20,7 +20,7 @@ set(SOURCE_FILES
         ./src/HomogeneousForcing.cpp ./src/HomogeneousForcing.h
         ./src/FunctionArguments.cpp ./src/FunctionArguments.h
         ./src/FiniteElementMeshInterface.cpp ./src/FiniteElementMeshInterface.h
-        )
+        src/SpatialForcing.cpp src/SpatialForcing.h)
 
 add_library(physics SHARED ${SOURCE_FILES})
 

src/Physics/src/FiniteElementMesh.cpp

@@ -78,6 +78,13 @@ std::vector<std::vector<XY>> FiniteElementMesh<P, Q>::quadrature_points() const
     return qp;
 }
 
+template<size_t P, size_t Q>
+std::vector<XY> FiniteElementMesh<P, Q>::quadrature_points(size_t i) const {
+    std::vector<XY> qp = Triangles[i].template quadrature_points<Q>(Nodes);
+
+    return qp;
+}
+
 // TODO: std::shared_ptr<DiscreteRegion> region(std::shared_ptr<Region const> const &r) const override {...}
 
 template<size_t P, size_t Q>

src/Physics/src/FiniteElementMesh.h

@@ -36,6 +36,8 @@ public:
 
     std::vector<std::vector<XY>> quadrature_points() const override;
 
+    std::vector<XY> quadrature_points(size_t i) const override;
+
     // TODO: std::shared_ptr<DiscreteRegion> region(std::shared_ptr<Region const> const &r) const override {...}
 
     std::vector<std::shared_ptr<DiscreteBoundary>> make_discontinuous(std::shared_ptr<Curve const> const &c) override;

src/Physics/src/FiniteElementMeshInterface.h

@@ -55,6 +55,8 @@ public:
 
     virtual std::vector<std::vector<XY>> quadrature_points() const = 0;
 
+    virtual std::vector<XY> quadrature_points(size_t i) const = 0;
+
     size_t nodes_size() const { return Nodes.size(); };
 
     virtual size_t elements_size() const = 0;

src/Physics/src/Forcing.h

@@ -5,12 +5,13 @@
 
 #include "DiscreteRegion.h"
 #include "FunctionArguments.h"
+#include "FiniteElementMeshInterface.h"
 
 class Forcing {
 public:
     virtual ~Forcing() {};
 
-    virtual Oe::VectorXd operator()(FunctionArguments fargs) const = 0;
+    virtual Oe::VectorXd operator()(FunctionArguments fargs, std::shared_ptr<FiniteElementMeshInterface> fem) const = 0;
 };
 
 #endif //OERSTED_CONDITION_H

src/Physics/src/HomogeneousForcing.cpp

 #include "HomogeneousForcing.h"
 
-HomogeneousForcing::HomogeneousForcing(std::function<double(FunctionArguments)> f, std::vector<std::shared_ptr<DiscreteRegion>> r, SparseMatrixGroup const &b,
-                                       DiagonalMatrixGroup const &w) : Function{f}, Regions{r}, Basis(b), Weights(w) {
+HomogeneousForcing::HomogeneousForcing(std::function<double_t(FunctionArguments)> f, std::vector<std::shared_ptr<DiscreteRegion>> r, SparseMatrixGroup const &b, DiagonalMatrixGroup const &w) : Function{f}, Regions{r}, Basis(b), Weights(w) {
     // TODO: May need to adjust implementation once quadrilaterals are added
     Indicator.resize(Weights.size());
     Indicator.setZero();
@@ -12,7 +11,7 @@ HomogeneousForcing::HomogeneousForcing(std::function<double(FunctionArguments)>
     }
 }
 
-Oe::VectorXd HomogeneousForcing::operator()(FunctionArguments fargs) const {
+Oe::VectorXd HomogeneousForcing::operator()(FunctionArguments fargs, std::shared_ptr<FiniteElementMeshInterface> fem) const {
     Oe::VectorXd output = (Basis(0) * (Weights(0).matrix().asDiagonal() * Indicator));
 
     for (size_t i = 1; i != Basis.size(); ++i) {

src/Physics/src/HomogeneousForcing.h

@@ -5,15 +5,14 @@
 
 class HomogeneousForcing : public Forcing {
 public:
-    HomogeneousForcing(std::function<double(FunctionArguments)> f, std::vector<std::shared_ptr<DiscreteRegion>> regions, SparseMatrixGroup const &basis,
-                       DiagonalMatrixGroup const &weights);
+    HomogeneousForcing(std::function<double_t(FunctionArguments)> f, std::vector<std::shared_ptr<DiscreteRegion>> regions, SparseMatrixGroup const &basis, DiagonalMatrixGroup const &weights);
 
-    Oe::VectorXd operator()(FunctionArguments fargs) const override;
+    Oe::VectorXd operator()(FunctionArguments fargs, std::shared_ptr<FiniteElementMeshInterface> fem) const override;
 
     auto const &indicator() const { return Indicator; };
 
 protected:
-    std::function<double(FunctionArguments)> Function;
+    std::function<double_t(FunctionArguments)> Function;
 
     std::vector<std::shared_ptr<DiscreteRegion>> Regions;
 

src/Physics/src/Physics.h

@@ -14,6 +14,7 @@
 #include "PostProcessorInterface.h"
 #include "SlidingInterface.h"
 #include "ZeroDirichlet.h"
+#include "SpatialForcing.h"
 
 // TODO: These classes should be generic (e.g. CurlCurl, ScalarLaplacian, VectorLaplacian)
 // TODO:    Higher level classes should implement interfaces with appropriate physics grammar
@@ -192,13 +193,13 @@ public:
     void calculate_forcing(Oe::VectorXd &f, FunctionArguments fargs) override {
         f.setZero();
         for (size_t i = 0; i != this->ForcingCondtions.size(); ++i) {
-            f += this->ForcingCondtions[i]->operator()(fargs);
+            f += this->ForcingCondtions[i]->operator()(fargs, Domain);
         }
     }
 
     // TODO: Make generic "add" methods in base class and provide interface with proper physics grammar in sub-classes (this can use forwarding to non-template base class impls)
 
-    void add_current_density(std::function<double(FunctionArguments)> func, std::vector<std::shared_ptr<Contour const>> contours) {
+    void add_current_density(std::function<double_t(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));
@@ -207,10 +208,23 @@ public:
         this->ForcingCondtions.push_back(std::make_shared<HomogeneousForcing>(func, regions, this->Basis, this->ElementWeights));
     }
 
-    void add_current_density(std::function<double(FunctionArguments)> func, std::vector<std::shared_ptr<DiscreteRegion>> regions) {
+    void add_current_density(std::function<double_t(FunctionArguments)> func, std::vector<std::shared_ptr<DiscreteRegion>> regions) {
         this->ForcingCondtions.push_back(std::make_shared<HomogeneousForcing>(func, regions, this->Basis, this->ElementWeights));
     }
 
+    void add_current_density(std::function<double_t(FunctionArguments,XY)> 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));
+        }
+
+        this->ForcingCondtions.push_back(std::make_shared<SpatialForcing>(func, regions, this->Basis, this->ElementWeights));
+    }
+
+    void add_current_density(std::function<double_t(FunctionArguments,XY)> func, std::vector<std::shared_ptr<DiscreteRegion>> regions) {
+        this->ForcingCondtions.push_back(std::make_shared<SpatialForcing>(func, regions, this->Basis, this->ElementWeights));
+    }
+
     void remove_current_density(size_t i) {
         this->ForcingCondtions.erase(this->ForcingCondtions.begin() + i);
     }

src/Physics/src/SpatialForcing.cpp

+#include "SpatialForcing.h"
+
+SpatialForcing::SpatialForcing(std::function<double_t(FunctionArguments,XY)> f, std::vector<std::shared_ptr<DiscreteRegion>> r, SparseMatrixGroup const &b, DiagonalMatrixGroup const &w) : Function{f}, Regions{r}, Basis(b), Weights(w) {
+    // TODO: May need to adjust implementation once quadrilaterals are added
+    for (auto const &r : Regions) {
+        for (auto const &i : r->elements()) {
+            Elements.push_back(i);
+        }
+    }
+}
+
+Oe::VectorXd SpatialForcing::operator()(FunctionArguments fargs, std::shared_ptr<FiniteElementMeshInterface> fem) const {
+    // TODO: Make this default implementation and then specialize for simpler cases?
+    size_t const Q = Basis.size();
+    size_t const N = (size_t)Weights(0).size();
+
+    Oe::MatrixXd qvals  = Oe::MatrixXd::Zero(Q, N);
+    for (auto const &n : Elements) {
+        auto qpts = fem->quadrature_points(n);
+        for (size_t q = 0; q != Q; ++q) {
+            qvals(q,n) += Function(fargs,qpts[q]);
+        }
+    }
+    qvals.transposeInPlace();
+
+    Oe::VectorXd output = Oe::VectorXd::Zero(Basis(0).rows());
+    for (size_t q = 0; q != Q; ++q) {
+        output += (Basis(q) * (Weights(q).matrix().asDiagonal() * qvals.col(q)));
+    }
+
+    return output;
+}

src/Physics/src/SpatialForcing.h

+#ifndef OERSTED_SPATIALFORCING_H
+#define OERSTED_SPATIALFORCING_H
+
+#include "Elements.hpp"
+
+#include "Forcing.h"
+
+class SpatialForcing : public Forcing {
+public:
+    SpatialForcing(std::function<double_t(FunctionArguments,XY)> f, std::vector<std::shared_ptr<DiscreteRegion>> regions, SparseMatrixGroup const &basis, DiagonalMatrixGroup const &weights);
+
+    Oe::VectorXd operator()(FunctionArguments fargs, std::shared_ptr<FiniteElementMeshInterface> fem) const override;
+
+protected:
+    std::function<double_t(FunctionArguments,XY)> Function;
+
+    std::vector<std::shared_ptr<DiscreteRegion>> Regions;
+
+    SparseMatrixGroup const &Basis;
+    DiagonalMatrixGroup const &Weights;
+
+    std::vector<size_t> Elements; // TODO?: Make Indicator part of the Region class instead? ElementIndicator, NodeIndicator
+};
+
+#endif //OERSTED_SPATIALFORCING_H

test/CMakeLists.txt

@@ -34,11 +34,11 @@ set(SOURCE_FILES
         LibraryIntegration/Mesh_To_FEM_Test.cpp
         LibraryIntegration/test_Square.cpp
         LibraryIntegration/test_Circle.cpp
+        LibraryIntegration/test_Annuli.cpp
 
         UseCases/Use_Cases.hpp
         UseCases/Stator_Use_Cases.cpp
-        UseCases/Rotor_Use_Cases.cpp
-        )
+        UseCases/Rotor_Use_Cases.cpp)
 
 add_executable(run_tests ${SOURCE_FILES})
 

test/LibraryIntegration/test_Annuli.cpp

+#include "Library_Integration_Test.hpp"
+
+class Annuli : public ::testing::Test {
+public:
+    virtual void SetUp() {
+        SetUp_Sketch();
+        SetUp_Mesh();
+    }
+
+    void SetUp_Sketch() {
+        auto v0 = sk.new_element<Vertex>(0.0, 0.0);
+
+        auto v1 = sk.new_element<Vertex>(+1.0, 0.0);
+        auto v2 = sk.new_element<Vertex>(0.0, +1.0);
+        auto v3 = sk.new_element<Vertex>(-1.0, 0.0);
+        auto v4 = sk.new_element<Vertex>(0.0, -1.0);
+
+        auto v5 = sk.new_element<Vertex>(+2.0, 0.0);
+        auto v6 = sk.new_element<Vertex>(0.0, +2.0);
+        auto v7 = sk.new_element<Vertex>(-2.0, 0.0);
+        auto v8 = sk.new_element<Vertex>(0.0, -2.0);
+
+        sk.new_element<CircularArc>(v1,v2,v0,1.0);
+        sk.new_element<CircularArc>(v2,v3,v0,1.0);
+        sk.new_element<CircularArc>(v3,v4,v0,1.0);
+        sk.new_element<CircularArc>(v4,v1,v0,1.0);
+
+        boundaries.push_back(sk.new_element<CircularArc>(v5,v6,v0,2.0));
+        boundaries.push_back(sk.new_element<CircularArc>(v6,v7,v0,2.0));
+        boundaries.push_back(sk.new_element<CircularArc>(v7,v8,v0,2.0));
+        boundaries.push_back(sk.new_element<CircularArc>(v8,v5,v0,2.0));
+
+        double residual = sk.solve();
+        EXPECT_LE(residual, FLT_EPSILON);
+
+        bool result = sk.build();
+        ASSERT_TRUE(result);
+
+        sk.save_as<SaveMethod::Rasterize>(SAVE_DIR, std::string("sketch"));
+    }
+
+    void SetUp_Mesh() {
+        me = Mesh(sk);
+
+        me.create();
+
+        me.MaximumElementSize = 0.1;
+        me.MinimumElementSize = 0.01;
+        me.MinimumElementQuality = 0.5;
+
+        bool result = me.refine();
+        ASSERT_TRUE(result);
+
+        me.save_as(SAVE_DIR, std::string("mesh"));
+
+        inner_region = me.select_contour(Point{0.0,0.0});
+        outer_region = me.select_contour(Point{1.5,0.0});
+    }
+
+    std::string SAVE_DIR = "./test/output/Integration/test_Annuli/";
+
+    std::vector<std::shared_ptr<Curve>> boundaries;
+
+    std::shared_ptr<const Contour> inner_region;
+    std::shared_ptr<const Contour> outer_region;
+
+    Sketch sk;
+    Mesh me;
+
+    std::shared_ptr<FiniteElementMeshInterface> fem;
+
+    // Homogenous Inner Forcing
+    template<size_t P, size_t Q>
+    auto magnetostatic_inner_homogenous_forcing() {
+        fem = std::make_shared<FiniteElementMesh<P, Q>>(me);
+        Magnetostatic ms{fem};
+
+        FunctionArguments fargs;
+
+        ms.add_current_density([](FunctionArguments args) { return 1.0; }, {inner_region});
+
+        ms.add_magnetic_insulation({boundaries[0],boundaries[1],boundaries[2],boundaries[3]});
+
+        ms.assemble();
+
+        auto sol = ms.new_solution();
+
+        ms.solve(sol, fargs);
+
+        auto v = ms.recover_boundary(sol->v);
+
+        return v;
+    };
+
+    void test_mihf_solution(Oe::VectorXd const &v, double_t tol) {
+        // v(r) = -mu_0 * (0.25 * r^2 + 0.5 * log(0.5) - 0.25), r <= 1.0
+        // v(r) = -mu_0 * 0.5 * log(0.5 * r),1.0 <= r <= 2.0
+
+        tol *= mu_0 * (0.25 - 0.5 * log(0.5));
+        for (size_t i = 0; i != fem->nodes().size(); ++i) {
+            double_t r = std::hypot(fem->node(i).x(),fem->node(i).y());
+
+            if (r < 1.0) {
+                double_t v_expected = -mu_0 * (0.25 * r * r + 0.5 * log(0.5) - 0.25);
+                EXPECT_NEAR(v(i), v_expected, tol);
+            } else {
+                double_t v_expected = -mu_0 * 0.5 * log(0.5 * r);
+                EXPECT_NEAR(v(i), v_expected, tol);
+            }
+        }
+    }
+
+    // Homogenous Outer Forcing
+    template<size_t P, size_t Q>
+    auto magnetostatic_outer_homogenous_forcing() {
+        fem = std::make_shared<FiniteElementMesh<P, Q>>(me);
+        Magnetostatic ms{fem};
+
+        FunctionArguments fargs;
+
+        ms.add_current_density([](FunctionArguments args) { return 1.0; }, {outer_region});
+
+        ms.add_magnetic_insulation({boundaries[0],boundaries[1],boundaries[2],boundaries[3]});
+
+        ms.assemble();
+
+        auto sol = ms.new_solution();
+
+        ms.solve(sol, fargs);
+
+        auto v = ms.recover_boundary(sol->v);
+
+        return v;
+    };
+
+    void test_mohf_solution(Oe::VectorXd const &v, double_t tol) {
+        // v(r) = -mu_0 * (-0.75 - 0.5 * log(0.5)), r <= 1.0
+        // v(r) = -mu_0 * (0.25 * r^2 - 0.5 * log(0.5 * r) - 1) ,1.0 <= r <= 2.0
+
+        tol *= mu_0 * (0.75 + 0.5 * log(0.5));
+        for (size_t i = 0; i != fem->nodes().size(); ++i) {
+            double_t r = std::hypot(fem->node(i).x(),fem->node(i).y());
+
+            if (r < 1.0) {
+                double_t v_expected = -mu_0 * (-0.75 - 0.5 * log(0.5));
+                EXPECT_NEAR(v(i), v_expected, tol);
+            } else {
+                double_t v_expected = -mu_0 * (0.25 * r * r - 0.5 * log(0.5 * r) - 1.0);
+                EXPECT_NEAR(v(i), v_expected, tol);
+            }
+        }
+    }
+
+    // Sinusoidal Inner Forcing
+    template<size_t P, size_t Q>
+    auto magnetostatic_inner_sinusoidal_forcing() {
+        fem = std::make_shared<FiniteElementMesh<P, Q>>(me);
+        Magnetostatic ms{fem};
+
+        FunctionArguments fargs;
+        auto func = [](FunctionArguments args, XY p) {
+            return std::sin(std::atan2(p.y(), p.x()));
+        };
+
+        ms.add_current_density(func, {inner_region});
+
+        ms.add_magnetic_insulation({boundaries[0],boundaries[1],boundaries[2],boundaries[3]});
+
+        ms.assemble();
+
+        auto sol = ms.new_solution();
+
+        ms.solve(sol, fargs);
+
+        auto v = ms.recover_boundary(sol->v);
+
+        return v;
+    };
+
+    void test_misf_solution(Oe::VectorXd const &v, double_t tol) {
+        // v(r) = -mu_0 * sin(a) * (r^2 / 3 - 11 * r / 24), r <= 1.0
+        // v(r) = -mu_0 * sin(a) * (r / 24 - 1 / (6 * r)) ,1.0 <= r <= 2.0
+
+        tol *= mu_0 * (1.0 / 8.0);
+        for (size_t i = 0; i != fem->nodes().size(); ++i) {
+            double_t r = std::hypot(fem->node(i).x(),fem->node(i).y());
+            double_t a = std::atan2(fem->node(i).y(),fem->node(i).x());
+
+            if (r < 1.0) {
+                double_t v_expected = -mu_0 * std::sin(a) * (r * r / 3.0 - 11.0 * r / 24.0);
+                EXPECT_NEAR(v(i), v_expected, tol);
+            } else {
+                double_t v_expected = -mu_0 * std::sin(a) * (r / 24.0 - 1.0 / (6.0 * r));
+                EXPECT_NEAR(v(i), v_expected, tol);
+            }
+        }
+    }
+};
+
+TEST_F(Annuli, magnetostatic_inner_homogenous_forcing) {
+    test_mihf_solution(magnetostatic_inner_homogenous_forcing<1,1>(), 0.01);
+    test_mihf_solution(magnetostatic_inner_homogenous_forcing<1,2>(), 0.01);
+
+    test_mihf_solution(magnetostatic_inner_homogenous_forcing<2,2>(), 0.0001); // Does not obtain FLT_EPSILON accuracy because of boundary discretization error
+}
+
+TEST_F(Annuli, magnetostatic_outer_homogenous_forcing) {
+    test_mohf_solution(magnetostatic_outer_homogenous_forcing<1,1>(), 0.01);
+    test_mohf_solution(magnetostatic_outer_homogenous_forcing<1,2>(), 0.01);
+
+    test_mohf_solution(magnetostatic_outer_homogenous_forcing<2,2>(), 0.0001); // Does not obtain FLT_EPSILON accuracy because of boundary discretization error
+}
+
+TEST_F(Annuli, magnetostatic_inner_sinusoidal_forcing) {
+    test_misf_solution(magnetostatic_inner_sinusoidal_forcing<1,1>(), 0.1);
+    test_misf_solution(magnetostatic_inner_sinusoidal_forcing<1,2>(), 0.01);
+
+
+    test_misf_solution(magnetostatic_inner_sinusoidal_forcing<2,2>(), 0.01);
+}
+
+TEST_F(Annuli, magnetostatic_outer_sinusoidal_forcing) {
+    FAIL(); // TODO
+}
+
+TEST_F(Annuli, magnetostatic_combined_sinusoidal_forcing) {
+    FAIL(); // TODO
+}
+
+TEST_F(Annuli, TODO) {
+    FAIL();
+
+    // TODO: Implement homogenous and sinusoidal tests on two coupled annuli (test_Annuli)
+    // TODO: Implement torque calculation using Maxwell stress tensor
+    // TODO: Test torque calculation in test_Square, test_Circle (both should be zero torque), and test_Annuli (certain configurations should give non-zero torque)
+
+    // TODO: Fix implementation so tests do not sigbart when "run all" is selected"
+}
\ No newline at end of file

test/LibraryIntegration/test_Circle.cpp

@@ -56,8 +56,9 @@ public:
 
     std::shared_ptr<FiniteElementMeshInterface> fem;
 
+    // Homogenous Forcing
     template<size_t P, size_t Q>
-    auto magnetostatic() {
+    auto magnetostatic_homogenous_forcing() {
         fem = std::make_shared<FiniteElementMesh<P, Q>>(me);
         Magnetostatic ms{fem};
 
@@ -78,7 +79,7 @@ public:
         return v;
     };
 
-    void test_solution(Oe::VectorXd const &v, double_t tol) {
+    void test_mhf_solution(Oe::VectorXd const &v, double_t tol) {
         // v(r) = a * (1 - r^2), a = mu_0 / 4.0;
 
         double_t a = mu_0 / 4.0;
@@ -88,17 +89,58 @@ public:
             EXPECT_NEAR(v(i), v_expected, tol * a);
         }
     }
+
+    // Sinusoidal Forcing
+    template<size_t P, size_t Q>
+    auto magnetostatic_sinusoidal_forcing() {
+        fem = std::make_shared<FiniteElementMesh<P, Q>>(me);
+        Magnetostatic ms{fem};
+
+        FunctionArguments fargs;
+        auto func = [](FunctionArguments args, XY p) {
+            return std::sin(std::atan2(p.y(), p.x()));
+        };
+
+        ms.add_current_density(func, {circle_contour});
+
+        ms.add_magnetic_insulation({boundaries[0],boundaries[1],boundaries[2],boundaries[3]});
+
+        ms.assemble();
+
+        auto sol = ms.new_solution();
+
+        ms.solve(sol, fargs);
+
+        auto v = ms.recover_boundary(sol->v);
+
+        return v;
+    };
+
+    void test_msf_solution(Oe::VectorXd const &v, double_t tol) {
+        // v(r,a) = mu_0 * (r - r^2) * sin(a) / 3
+        // max(v(r,a)) = mu_0 / 12.0 @ (r,a) = (1/2,pi/2)
+
+        for (size_t i = 0; i != fem->nodes().size(); ++i) {
+            double_t r = std::hypot(fem->node(i).x(),fem->node(i).y());
+            double_t a = std::atan2(fem->node(i).y(),fem->node(i).x());
+
+            double_t v_expected = mu_0 * (r - r * r) * sin(a) / 3.0;
+
+            EXPECT_NEAR(v(i), v_expected, tol * mu_0 / 12.0);
+        }
+    }
 };
 
-TEST_F(Circle, magnetostatic) {
-    test_solution(magnetostatic<1,1>(), 0.01);
-    test_solution(magnetostatic<1,2>(), 0.01);
-    test_solution(magnetostatic<2,2>(), 0.0001); // Does not obtain FLT_EPSILON accuracy because of boundary discretization error
-
-    FAIL();
-    // TODO: Implement spatially dependent current distribution forcing function
-    // TODO: Implement tests on circle with sinusoidal current distribution (cylindrical coordinates)
-    // TODO: Implement tests on two coupled annuli (test_Annuli)
-    // TODO: Implement torque calculation using Maxwell stress tensor
-    // TODO: Test torque calculation in test_Square, test_Circle (both should be zero torque), and test_Annuli (certain configurations should give non-zero torque)
+TEST_F(Circle, magnetostatic_homogenous_forcing) {
+    test_mhf_solution(magnetostatic_homogenous_forcing<1,1>(), 0.01);
+    test_mhf_solution(magnetostatic_homogenous_forcing<1,2>(), 0.01);
+
+    test_mhf_solution(magnetostatic_homogenous_forcing<2,2>(), 0.0001); // Does not obtain FLT_EPSILON accuracy because of boundary discretization error
+}
+
+TEST_F(Circle, magnetostatic_sinusoidal_forcing) {
+    test_msf_solution(magnetostatic_sinusoidal_forcing<1,1>(), 0.1);
+    test_msf_solution(magnetostatic_sinusoidal_forcing<1,2>(), 0.1);
+
+    test_msf_solution(magnetostatic_sinusoidal_forcing<2,2>(), 0.01);
 }
\ No newline at end of file

test/LibraryIntegration/test_Square.cpp

@@ -56,8 +56,9 @@ public:
 
     std::shared_ptr<FiniteElementMeshInterface> fem;
 
+    // Homogenous Forcing
     template<size_t P, size_t Q>
-    auto magnetostatic(size_t b) {
+    auto magnetostatic_homogenous_forcing(size_t b) {
         fem = std::make_shared<FiniteElementMesh<P, Q>>(me);
         Magnetostatic ms{fem};
 
@@ -78,7 +79,7 @@ public:
         return v;
     };
 
-    void test_boundary0_solution(Oe::VectorXd const &v, double_t tol) {
+    void test_mhf_boundary0_solution(Oe::VectorXd const &v, double_t tol) {
         // v(x,y) = -0.5 * mu_o * y^2 + mu_o * y
         for (size_t i = 0; i != fem->nodes().size(); ++i) {
             double_t y = fem->node(i).y();
@@ -87,7 +88,7 @@ public:
         }
     }
 
-    void test_boundary3_solution(Oe::VectorXd const &v, double_t tol) {
+    void test_mhf_boundary3_solution(Oe::VectorXd const &v, double_t tol) {
         // v(x,y) = -0.5 * mu_o * x^2 + mu_o * x
         for (size_t i = 0; i != fem->nodes().size(); ++i) {
             double_t x = fem->node(i).x();
@@ -95,14 +96,61 @@ public:
             EXPECT_NEAR(v(i), v_expected, tol * 0.5 * mu_0);
         }
     }
+