Mesh::mesh_data() now returns a MeshData object containing a vector of...

Mesh::mesh_data() now returns a MeshData object containing a vector of Triangles instead of nested vectors
Added SquareOrder2 tests of second order triangular elements

src/Elements/CMakeLists.txt

@@ -5,9 +5,9 @@ set(SOURCE_FILES
 
         ./src/Triangle.cpp  ./src/Triangle.h
         ./src/Node.cpp      ./src/Node.h
+        ./src/Point.cpp     ./src/Point.h
         ./src/Element.cpp   ./src/Element.h
-        ./src/Facet.cpp     ./src/Facet.h
-        ../Materials/include/Materials.hpp)
+        ./src/Facet.cpp     ./src/Facet.h)
 
 add_library(elements SHARED ${SOURCE_FILES})
 

src/Elements/include/Elements.hpp

@@ -2,6 +2,7 @@
 #define OERSTED_ELEMENTS_HP
 
 #include "../src/Node.h"
+#include "../src/Point.h"
 #include "../src/Triangle.h"
 
 #endif //OERSTED_ELEMENTS_HPP

src/Elements/src/Node.h

@@ -4,10 +4,9 @@
 // TODO: Const
 
 #include <cstddef>
+#include <cmath>
 
-#include "Mesh.hpp"
-
-// TODO: Combine different node/point classes from Mesh and Elements
+#include "Point.h" // TODO: Combine different node/point classes from Mesh and Elements
 
 class XY {
 public:

src/Mesh/src/Point.cpp → src/Elements/src/Point.cpp

-#include "Mesh.hpp"
+#include "Point.h"
 
 double dist(Point const &p0, Point const &p1) {
     return std::sqrt((p0.X - p1.X) * (p0.X - p1.X) + (p0.Y - p1.Y) * (p0.Y - p1.Y));

src/Mesh/CMakeLists.txt

@@ -9,7 +9,6 @@ set(SOURCE_FILES
         ./src/InsertionQueuer.h     ./src/InsertionQueuer.cpp
         ./src/MappedBoundaryPair.h  ./src/MappedBoundaryPair.cpp
         ./src/Mesh.h                ./src/Mesh.cpp
-        ./src/Point.h               ./src/Point.cpp
         ./src/DartTriangle.h        ./src/DartTriangle.cpp)
 
 add_library(mesh SHARED ${SOURCE_FILES})

src/Mesh/src/Edge.h

 #ifndef OERSTED_EDGE_H
 #define OERSTED_EDGE_H
 
+#include "Elements.hpp"
+
 #include "DartConstraint.h"
-#include "Point.h"
 #include "InsertionQueuer.h"
 
 class Mesh;

src/Mesh/src/InsertionQueuer.h

@@ -3,7 +3,7 @@
 
 #include <cstddef>
 
-#include "Point.h"
+#include "Elements.hpp"
 
 class Mesh;
 

src/Mesh/src/Mesh.h

@@ -11,12 +11,12 @@
 #include "Eigen/SparseLU"
 
 #include "Sketch.hpp"
+#include "Elements.hpp"
 
 #include "BoundaryConstraint.h"
 #include "MappedBoundaryPair.h"
 #include "Edge.h"
 #include "InsertionQueuer.h"
-#include "Point.h"
 #include "DartTriangle.h"
 
 enum class LocateTriangleResult {
@@ -48,11 +48,14 @@ public:
     Point Coordinates;
 };
 
-class MeshData {
+template<size_t P>
+class MeshData{
 public:
     std::vector<Point> Nodes;
 
-    std::vector<std::vector<std::vector<size_t>>> Edges;
+    // std::vector<std::vector<std::vector<size_t>>> Edges;
+
+    std::vector<Triangle<P>> Triangles;
 
     std::map<std::shared_ptr<const Curve>, std::set<std::pair<double_t,size_t>>> Boundaries;
 
@@ -161,8 +164,7 @@ public:
         /*
          * Constructs a MappedBoundaryPair mesh constraint from a ContinuousBoundaryPair
          *
-         * //TODO: Present implementation strategy allows cyclical mapping
-         * //TODO: which should converge for boundaries restricted to power-of-2 parametric discretizations, but is untested
+         * //TODO: Present implementation strategy allows cyclical mapping which should converge for boundaries restricted to power-of-2 parametric discretizations, but is untested
          */
 
         MappedBoundaries.push_back(std::make_shared<MappedBoundaryPair>(*this, cbp));
@@ -233,12 +235,11 @@ public:
         return n;
     };
 
-    MeshData mesh_data(size_t p) const {
-        // TODO: Modify this to return triangles and nodes instead of std::vector<std::vector<std::vector<size_t>>>
-
+    template<size_t P>
+    MeshData<P> mesh_data() const {
         std::vector<NodeData> data;
         std::vector<Point> points;
-        std::vector<std::vector<std::vector<size_t>>> edges;
+        std::vector<Triangle<P>> triangles;
 
         std::map<std::shared_ptr<const Curve>, std::set<std::pair<double_t,size_t>>> boundaries;
         std::map<std::shared_ptr<const Contour>, std::vector<size_t>> regions;
@@ -261,10 +262,10 @@ public:
                     param1 = dc.S0;
                 }
 
-                for (size_t i = 0; i != p + 1; i++) {
-                    double_t s = i / (double_t) (p);
-                    Point pointi = point0 * s + point1 * (1.0 - s);
-                    double_t parami = param0 * s + param1 * (1.0 - s);
+                for (size_t i = 0; i != P + 1; i++) {
+                    double_t s = i / (double_t) (P);
+                    Point pointi = point1 * s + point0 * (1.0 - s);
+                    double_t parami = param1 * s + param0 * (1.0 - s);
 
                     data.emplace_back(t, k, i, dc.curve(), parami, pointi);
                 }
@@ -272,28 +273,25 @@ public:
                 e = next(e);
             }
         }
-
         std::sort(data.begin(), data.end());
 
-        edges.resize(p + 1);
-        for (size_t i = 0; i != p + 1; ++i) {
-            edges[i].resize(3);
-            for (size_t j = 0; j != 3; ++j) {
-                edges[i][j].resize(size_triangles());
-            }
-        }
+        triangles.resize(size_triangles());
 
         size_t i = 0;
         size_t n = 0;
         while (i < data.size()) {
             points.push_back(data[i].Coordinates);
 
-            edges[data[i].NodeID][data[i].EdgeID][data[i].ElementID] = n;
+            size_t k = data[i].EdgeID * P + data[i].NodeID;
+            triangles[data[i].ElementID].node(k,n);
+
             boundaries[data[i].Constraint].emplace(data[i].Parameter,n);
 
             size_t j = i + 1;
             while (j < data.size() && data[i].Coordinates == data[j].Coordinates) {
-                edges[data[j].NodeID][data[j].EdgeID][data[j].ElementID] = n;
+                size_t k = data[j].EdgeID * P + data[j].NodeID;
+                triangles[data[j].ElementID].node(k,n);
+
                 boundaries[data[j].Constraint].emplace(data[j].Parameter,n);
 
                 ++j;
@@ -302,9 +300,12 @@ public:
             i = j;
             ++n;
         }
-        edges.erase(edges.end()-1);
 
-        return MeshData{points,edges,boundaries,regions};
+        for(size_t i = 0; i != triangles.size(); ++i) {
+            triangles[i].id(i);
+        }
+
+        return MeshData<P>{points,triangles,boundaries,regions};
     }
 
     std::vector<size_t> boundary_nodes(size_t i) const { return BoundaryConstraints[i]->nodes(*this); }

src/Physics/src/FiniteElementMesh.cpp

@@ -13,41 +13,24 @@ FiniteElementMesh<P, Q>::FiniteElementMesh(Mesh &m) {
         Perhaps a similar approach can be taken for implementing periodic boundary conditions (e.g. MappedBoundary holds two (or several arrays of) boundaries and enforces the proper invariant)
             e.g. make_mapped(...)
         Then PeriodicBoundaryCondition holds a pointer to the MappedBoundary instead of two Boundary objects
-
-        TODO: To make these easier to implement from the user side, Boundary should hold a shared_ptr to the defining curve (from the Mesh object)
-            Then selection can be implemented by passing the shared_ptr that the user has (similar to how the Mesh operations work)
      */
 
-    // TODO: CLEAN ME and MESH
-
-    MeshData md = m.mesh_data(P);
+    MeshData<P> md = m.mesh_data<P>();
 
     Nodes.reserve(md.Nodes.size());
     for (size_t i = 0; i != md.Nodes.size(); ++i) {
         Nodes.emplace_back(md.Nodes[i]);
     }
 
-    //Boundaries.reserve(md.Boundaries.size());
     for (auto b = ++md.Boundaries.begin(); b != md.Boundaries.end(); ++b) {
         Boundaries.emplace(b->first, std::make_shared<DiscreteBoundary>(b->first, b->second));
     }
 
-    //Regions.reserve(md.Regions.size());
     for (auto r = md.Regions.begin(); r != md.Regions.end(); ++r) {
         Regions.emplace(r->first, std::make_shared<DiscreteRegion>(r->first, r->second)); // TODO: Assign material properties here
     }
 
-    Triangles.resize(md.Edges[0][0].size());
-    size_t n = 0;
-    for (size_t i = 0; i != md.Edges.size(); ++i) {
-        for (size_t j = 0; j != md.Edges[i].size(); ++j) {
-            for (size_t k = 0; k != md.Edges[i][j].size(); ++k) {
-                Triangles[k].id(k); //TODO: FIX THIS STUPID THING
-                Triangles[k].node(n, md.Edges[i][j][k]);
-            }
-            ++n;
-        }
-    };
+    Triangles = md.Triangles;
 };
 
 template<size_t P, size_t Q>
@@ -105,7 +88,7 @@ std::vector<std::shared_ptr<DiscreteBoundary>> FiniteElementMesh<P, Q>::make_dis
     if (b_vec.size() == 0 || b_vec.size() == 2) { // Boundary curve not found || boundary pair found
         return b_vec;
     } else if (b_vec.size() > 2) { // should not happen
-        throw std::runtime_error("FiniteElementMesh<2>::make_discontinuous expects to find at most 2 DiscreteBoundary objects");
+        throw std::runtime_error("FiniteElementMesh::make_discontinuous expects to find at most 2 DiscreteBoundary objects");
     } // else b.size() == 1 and we need to construct an new boundary
 
     // Copy nodes
@@ -335,4 +318,7 @@ void FiniteElementMesh<P, Q>::write_vector(Oe::ArrayXd const &Fx, Oe::ArrayXd co
 
 // Explicit Instantiations
 template class FiniteElementMesh<1, 1>;
-template class FiniteElementMesh<1, 2>;
\ No newline at end of file
+template class FiniteElementMesh<1, 2>;
+
+template class FiniteElementMesh<2, 1>;
+template class FiniteElementMesh<2, 2>;
\ No newline at end of file

test/CMakeLists.txt

@@ -23,18 +23,23 @@ set(SOURCE_FILES
         Mesh/Mest_Test_Utility.h
         Mesh/Mest_Test_Utility.cpp
 
+        Elements/test_Elements.hpp
         Elements/test_Triangle.cpp
 
         Physics/test_Physics.hpp
         Physics/Finite_Element_Mesh_Test.cpp
+        Physics/test_SquareOrder2.cpp
 
         LibraryIntegration/Library_Integration_Test.hpp
         LibraryIntegration/Mesh_To_FEM_Test.cpp
+        LibraryIntegration/test_Square.cpp
+        LibraryIntegration/test_Circle.cpp
 
         UseCases/Use_Cases.hpp
         UseCases/Stator_Use_Cases.cpp
-        UseCases/Rotor_Use_Cases.cpp Elements/test_Triangle.cpp Elements/test_Elements.hpp)
+        UseCases/Rotor_Use_Cases.cpp
+        )
 
 add_executable(run_tests ${SOURCE_FILES})
 
-target_link_libraries(run_tests gtest gtest_main sketch mesh materials quadrature physics stdc++fs)
\ No newline at end of file
+target_link_libraries(run_tests gtest gtest_main sketch elements mesh materials quadrature physics stdc++fs)
\ No newline at end of file

test/LibraryIntegration/test_Circle.cpp

+#include "Library_Integration_Test.hpp"
+
+class Circle : 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>(1.0, 1.0);
+        auto v3 = sk.new_element<Vertex>(0.0, 1.0);
+
+        boundaries.push_back(sk.new_element<LineSegment>(v0,v1));
+        boundaries.push_back(sk.new_element<LineSegment>(v1,v2));
+        boundaries.push_back(sk.new_element<LineSegment>(v2,v3));
+        boundaries.push_back(sk.new_element<LineSegment>(v3,v0));
+
+        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() {
+        // Create Refineable 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"));
+
+        square_contour = me.select_contour(Point{0.5,0.5});
+    }
+
+    std::string SAVE_DIR = "./test/output/Integration/test_Square";
+
+    std::vector<std::shared_ptr<Curve>> boundaries;
+
+    std::shared_ptr<const Contour> square_contour;
+
+    Sketch sk;
+    Mesh me;
+
+    std::shared_ptr<FiniteElementMeshInterface> fem;
+
+    template<size_t P, size_t Q>
+    auto magnetostatic(size_t b) {
+        fem = std::make_shared<FiniteElementMesh<P, Q>>(me);
+        Magnetostatic ms{fem};
+
+        FunctionArguments fargs;
+
+        ms.add_current_density([](FunctionArguments args) { return 1.0; }, {square_contour});
+
+        ms.add_magnetic_insulation({boundaries[b]});
+
+        ms.assemble();
+
+        auto sol = ms.new_solution();
+
+        ms.solve(sol, fargs);
+
+        auto v = ms.recover_boundary(sol->v);
+
+        return v;
+    };
+
+    void test_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();
+            double_t v_expected = -0.5 * mu_0 * y * y + mu_0 * y;
+            EXPECT_NEAR(v(i), v_expected, tol * 0.5 * mu_0);
+        }
+    }
+
+    void test_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();
+            double_t v_expected = -0.5 * mu_0 * x * x + mu_0 * x;
+            EXPECT_NEAR(v(i), v_expected, tol * 0.5 * mu_0);
+        }
+    }
+    */
+};
+
+TEST_F(Circle, magnetostatic) {
+    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)
+}
\ No newline at end of file

test/LibraryIntegration/test_Square.cpp

+#include "Library_Integration_Test.hpp"
+
+class Square : 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>(1.0, 1.0);
+        auto v3 = sk.new_element<Vertex>(0.0, 1.0);
+
+        boundaries.push_back(sk.new_element<LineSegment>(v0,v1));
+        boundaries.push_back(sk.new_element<LineSegment>(v1,v2));
+        boundaries.push_back(sk.new_element<LineSegment>(v2,v3));
+        boundaries.push_back(sk.new_element<LineSegment>(v3,v0));
+
+        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() {
+        // Create Refineable 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"));
+
+        square_contour = me.select_contour(Point{0.5,0.5});
+    }
+
+    std::string SAVE_DIR = "./test/output/Integration/test_Square";
+
+    std::vector<std::shared_ptr<Curve>> boundaries;
+
+    std::shared_ptr<const Contour> square_contour;
+
+    Sketch sk;
+    Mesh me;
+
+    std::shared_ptr<FiniteElementMeshInterface> fem;
+
+    template<size_t P, size_t Q>
+    auto magnetostatic(size_t b) {
+        fem = std::make_shared<FiniteElementMesh<P, Q>>(me);
+        Magnetostatic ms{fem};
+
+        FunctionArguments fargs;
+
+        ms.add_current_density([](FunctionArguments args) { return 1.0; }, {square_contour});
+
+        ms.add_magnetic_insulation({boundaries[b]});
+
+        ms.assemble();
+
+        auto sol = ms.new_solution();
+
+        ms.solve(sol, fargs);
+
+        auto v = ms.recover_boundary(sol->v);
+
+        return v;
+    };
+
+    void test_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();
+            double_t v_expected = -0.5 * mu_0 * y * y + mu_0 * y;
+            EXPECT_NEAR(v(i), v_expected, tol * 0.5 * mu_0);
+        }
+    }
+
+    void test_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();
+            double_t v_expected = -0.5 * mu_0 * x * x + mu_0 * x;
+            EXPECT_NEAR(v(i), v_expected, tol * 0.5 * mu_0);
+        }
+    }
+};
+
+TEST_F(Square, magnetostatic) {
+    test_boundary0_solution(magnetostatic<1,1>(0), 0.01);
+    test_boundary0_solution(magnetostatic<1,2>(0), 0.01);
+    test_boundary0_solution(magnetostatic<2,2>(0), FLT_EPSILON);
+
+    test_boundary3_solution(magnetostatic<1,1>(3), 0.01);
+    test_boundary3_solution(magnetostatic<1,2>(3), 0.01);
+    test_boundary3_solution(magnetostatic<2,2>(3), FLT_EPSILON);
+}
\ No newline at end of file

test/Physics/Finite_Element_Mesh_Test.cpp

@@ -250,11 +250,11 @@ public:
         tri.push_back(Triangle<1>(0, {0, 1, 2}));
         tri.push_back(Triangle<1>(1, {1, 3, 2}));
 
-        reg.emplace(std::make_shared<Contour>(), std::make_shared<DiscreteRegion>(std::vector<size_t>{0}));
-        reg.emplace(std::make_shared<Contour>(), std::make_shared<DiscreteRegion>(std::vector<size_t>{1}));
-        for (auto it : reg) {
-            creg.push_back(it.first);
+        for (size_t i = 0; i != 2; ++i) {
+            creg.push_back(std::make_shared<Contour>());
         }
+        reg.emplace(creg[0], std::make_shared<DiscreteRegion>(std::vector<size_t>{0}));
+        reg.emplace(creg[1], std::make_shared<DiscreteRegion>(std::vector<size_t>{1}));
 
         bnd.emplace(std::make_shared<LineSegment>(), std::make_shared<DiscreteBoundary>(std::vector<size_t>{0, 1}));
         bnd.emplace(std::make_shared<LineSegment>(), std::make_shared<DiscreteBoundary>(std::vector<size_t>{0, 2}));

test/Physics/test_SquareOrder2.cpp

+#include "test_Physics.hpp"
+
+class SquareOrder2 : public ::testing::Test {
+public:
+    virtual void SetUp() {
+        node.push_back(XY{0.0, 0.0});//0
+        node.push_back(XY{0.5, 0.0});//1
+        node.push_back(XY{1.0, 0.0});//2
+        node.push_back(XY{0.0, 0.5});//3
+        node.push_back(XY{0.5, 0.5});//4
+        node.push_back(XY{1.0, 0.5});//5
+        node.push_back(XY{0.0, 1.0});//6
+        node.push_back(XY{0.5, 1.0});//7
+        node.push_back(XY{1.0, 1.0});//8
+
+        tri.push_back(Triangle<2>(0, {0,1,2,4,6,3}));
+        tri.push_back(Triangle<2>(1, {2,5,8,7,6,4}));
+
+        for (size_t i = 0; i != 2; ++i) {
+            creg.push_back(std::make_shared<Contour>());
+        }
+
+        reg.emplace(creg[0], std::make_shared<DiscreteRegion>(std::vector<size_t>{0}));
+        reg.emplace(creg[1], std::make_shared<DiscreteRegion>(std::vector<size_t>{1}));
+
+        for (size_t i = 0; i != 4; ++i) {
+            cbnd.push_back(std::make_shared<LineSegment>());
+        }
+
+        bnd.emplace(cbnd[0], std::make_shared<DiscreteBoundary>(std::vector<size_t>{0, 1, 2}));
+        bnd.emplace(cbnd[1], std::make_shared<DiscreteBoundary>(std::vector<size_t>{6, 7, 8}));
+        bnd.emplace(cbnd[2], std::make_shared<DiscreteBoundary>(std::vector<size_t>{0, 3, 6}));
+        bnd.emplace(cbnd[3], std::make_shared<DiscreteBoundary>(std::vector<size_t>{2, 5, 8}));
+
+        EXPECT_EQ(reg.size(), 2);
+        EXPECT_EQ(bnd.size(), 4);
+    }
+
+    std::vector<XY> node;
+    std::vector<Triangle<2>> tri;
+
+    std::map<std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> reg;
+    std::multimap<std::shared_ptr<const Curve>, std::shared_ptr<DiscreteBoundary>> bnd;
+
+    std::vector<std::shared_ptr<const Contour>> creg;
+    std::vector<std::shared_ptr<const Curve>> cbnd;
+
+    template<size_t QuadratureOrder>
+    std::shared_ptr<FiniteElementMesh<2, QuadratureOrder>> get_mesh() {
+        return std::make_shared<FiniteElementMesh<2,QuadratureOrder>>(node, tri, reg, bnd);
+    };
+
+    auto magnetostatic(size_t insulation_boundary_index) {
+        auto fem = get_mesh<2>();
+
+        std::shared_ptr<Magnetostatic> ms = std::make_shared<Magnetostatic>(fem);
+
+        // Add Current Density
+        FunctionArguments fargs;
+        fargs.add_argument("t", 1.0);
+        auto ff = [](FunctionArguments args) { return 1.0 * args["t"]; };
+        ms->add_current_density(ff, {reg[creg[0]], reg[creg[1]]});
+
+        // Add boundary condition
+        ms->add_magnetic_insulation({cbnd[insulation_boundary_index]});
+
+        // Solve
+        ms->assemble();
+        auto sol = ms->new_solution();
+
+        ms->solve(sol, fargs);
+
+        return sol;
+    }
+};
+
+TEST_F(SquareOrder2, magnetostatic) {
+    { // Boundary 0 is magnetic insulation
+        auto sol = magnetostatic(0);
+
+        // Test, v(x,y) = -0.5 * mu_o * y^2 + mu_o * y
+        std::vector<double_t> f_expected{1.0 / 6.0, 2.0 / 6.0, 1.0 / 6.0, 0.0, 1.0 / 6.0, 0.0};
+        std::vector<double_t> v_expected{3.0 / 8.0 * mu_0, 3.0 / 8.0 * mu_0, 3.0 / 8.0 * mu_0, 1.0 / 2.0 * mu_0, 1.0 / 2.0 * mu_0, 1.0 / 2.0 * mu_0};
+
+        for (size_t i = 0; i != 6; ++i) {
+            EXPECT_NEAR(so