Nightly Commit

Framework for multi-material magnetic models and testing of linear material properties with relative permeability > 1

src/Physics/CMakeLists.txt

@@ -9,6 +9,8 @@ set(SOURCE_FILES
 
         ./src/FiniteElementMesh.h ./src/FiniteElementMesh.cpp
 
+        ./src/MaterialProperties.h ./src/MaterialProperties.cpp
+
         ./src/MatrixGroup.h ./src/MatrixGroup.cpp
 
         ./src/Node.h ./src/Node.cpp

src/Physics/include/Physics.hpp

@@ -9,9 +9,11 @@
 #include "../src/BoundaryCondition.h"
 #include "../src/FiniteElementMesh.h"
 #include "../src/Forcing.h"
+#include "../src/MaterialProperties.h"
 #include "../src/MatrixGroup.h"
 #include "../src/Node.h"
 #include "../src/Physics.h"
+#include "../src/PhysicsConstants.h"
 #include "../src/QuadratureRule.h"
 #include "../src/Region.h"
 #include "../src/Triangle.h"

src/Physics/src/FiniteElementMesh.h

@@ -27,6 +27,7 @@ public:
 
     std::vector<Triangle<Order>> const &triangles() const { return Triangles; };
 
+    std::vector<Region<2>> &regions() { return Regions; };
     std::vector<Region<2>> const &regions() const { return Regions; };
 
     std::vector<Boundary<2>> const &boundaries() const { return Boundaries; };
@@ -35,6 +36,7 @@ public:
 
     Triangle<Order> const &triangle(size_t i) const { return Triangles[i]; };
 
+    Region<2> &region(size_t i) { return Regions[i]; };
     Region<2> const &region(size_t i) const { return Regions[i]; };
 
     Boundary<2> const &boundary(size_t i) const { return Boundaries[i]; };

src/Physics/src/Forcing.h

@@ -41,9 +41,9 @@ public:
     };
 
     Eigen::VectorXd operator()(double const t) const override {
-        Eigen::VectorXd output = (Basis[0] * (Weights(0) * Indicator)); // TODO: ?SparseVector?
+        Eigen::VectorXd output = (Basis[0] * (Weights(0).matrix().asDiagonal() * Indicator)); // TODO: ?SparseVector?
         for (size_t i = 1; i != TriangleQuadratureRule<QuadratureOrder>::size; ++i) {
-            output += (Basis[i] * (Weights(i) * Indicator));
+            output += (Basis[i] * (Weights(i).matrix().asDiagonal() * Indicator));
         }
 
         output *= Function(t);

src/Physics/src/MaterialProperties.cpp

+#include "MaterialProperties.h"

src/Physics/src/MaterialProperties.h

+#ifndef OERSTED_MATERIALPROPERTIES_H
+#define OERSTED_MATERIALPROPERTIES_H
+
+#include <memory>
+#include <utility>
+#include <vector>
+
+#include "Eigen"
+
+#include "PhysicsConstants.h"
+
+// TODO: Move into some sort of header
+using Vector = Eigen::VectorXd;
+using Array = Eigen::ArrayXd;
+
+class MagneticMaterialInterface {
+public:
+    virtual ~MagneticMaterialInterface(){};
+
+    virtual void h_from_b(std::vector<size_t> const &index, Array &Fx, Array &Fy, Array &dFxdx, Array &dFydy, Array &dFxdy) = 0;
+};
+
+class LinearIsotropicMagneticMaterial : public MagneticMaterialInterface {
+public:
+    LinearIsotropicMagneticMaterial(double relative_permeability) : Reluctivity(1.0 / (mu_0 * relative_permeability)) {};
+
+    void h_from_b(std::vector<size_t> const &index, Array &Fx, Array &Fy, Array &dFxdx, Array &dFydy, Array &dFxdy) override {
+        for(size_t const &i : index) {
+            Fx(i) *= Reluctivity;
+            Fy(i) *= Reluctivity;
+            dFxdx(i) = Reluctivity;
+            dFydy(i) = Reluctivity;
+            dFxdy(i) = 0.0;
+        }
+    }
+
+protected:
+    double Reluctivity;
+};
+
+class MaterialProperties {
+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} {};
+
+    template<typename... Args>
+    void h_from_b(Args&&... args) { return MagneticProperties->h_from_b(std::forward<Args>(args)...); }
+
+protected:
+    std::shared_ptr<MagneticMaterialInterface> MagneticProperties;
+};
+
+MaterialProperties Air() {
+    return MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1.0));
+}
+
+#endif //OERSTED_MATERIALPROPERTIES_H

src/Physics/src/MatrixGroup.h

@@ -47,7 +47,8 @@ public:
     auto &operator()(size_t const q) { return Matrices[q]; };
 
 protected:
-    std::array<Eigen::DiagonalMatrix<double, Eigen::Dynamic>, Q> Matrices;
+    //std::array<Eigen::DiagonalMatrix<double, Eigen::Dynamic>, Q> Matrices;
+    std::array<Eigen::ArrayXd, Q> Matrices;
 };
 
 template<size_t Q>

src/Physics/src/Physics.h

@@ -2,15 +2,18 @@
 #define OERSTED_PHYSICS_H
 
 #include <memory>
+#include <set>
 
-#include "Eigen/SparseCholesky"
 #include "Eigen/IterativeLinearSolvers"
+#include "Eigen/SparseCholesky"
 
-#include "Node.h"
-#include "Triangle.h"
-#include "QuadratureRule.h"
+#include "BoundaryCondition.h"
 #include "FiniteElementMesh.h"
 #include "Forcing.h"
+#include "Node.h"
+#include "PhysicsConstants.h"
+#include "QuadratureRule.h"
+#include "Triangle.h"
 
 enum class Variable {
     MagneticVectorPotential = 1,
@@ -102,37 +105,35 @@ public:
 
     Eigen::VectorXd init_element_vector() const { return Eigen::VectorXd(Elements); };
 
+    Eigen::ArrayXd init_unknown_array() const { return Eigen::ArrayXd(Unknowns); };
+
+    Eigen::ArrayXd init_element_array() const { return Eigen::ArrayXd(Elements); };
+
     void linearize(Eigen::SparseMatrix<double> &J, Eigen::VectorXd &v, Eigen::VectorXd &r, Eigen::VectorXd &f,
-                   Eigen::VectorXd &Fx, Eigen::VectorXd &Fy, Eigen::DiagonalMatrix<double, Eigen::Dynamic> &dFxdx, Eigen::DiagonalMatrix<double, Eigen::Dynamic> &dFydy, Eigen::DiagonalMatrix<double, Eigen::Dynamic> &dFxdy) {
+                   Eigen::ArrayXd &Fx, Eigen::ArrayXd &Fy, Eigen::ArrayXd &dFxdx, Eigen::ArrayXd &dFydy, Eigen::ArrayXd &dFxdy) {
 
         r = -f;
         J.setZero();
         for (size_t i = 0; i != TriangleQuadratureRule<QuadratureOrder>::size; ++i) {
             // Caclulate flux density
-            Fx = (Derivatives.dy(i).transpose() * v);
-            Fy = -(Derivatives.dx(i).transpose() * v);
+            Fx = Derivatives.dy(i).transpose() * v;
+            Fy = -Derivatives.dx(i).transpose() * v;
 
             // Calculate polarization
-            /*
-            for(auto const &r : domain.regions()) {
-                r.material().MagneticFieldIntensity(r.elements(),Fx,Fy,dFxdx,dFydy,dFxdy);
+            for(auto &dr : Domain.regions()) {
+                dr.material().h_from_b(dr.triangles(),Fx,Fy,dFxdx,dFydy,dFxdy);
             }
-            */
 
             // Calculate residual
-            r += (Derivatives.dy(i) * (ElementWeights(i) * Fx))
-                 - (Derivatives.dx(i) * (ElementWeights(i) * Fy)); // note: weak form introduces a negative sign into double-curl operator, hence r -= ...
+            r += Derivatives.dy(i) * (ElementWeights(i) * Fx).matrix()
+                 - Derivatives.dx(i) * (ElementWeights(i) * Fy).matrix(); // note: weak form introduces a negative sign into double-curl operator
 
             // Calculate jacobian
-            J += (Derivatives.dx(i) * ElementWeights(i) * Derivatives.dx(i).transpose())
-                 + (Derivatives.dy(i) * ElementWeights(i) * Derivatives.dy(i).transpose());
-
-            /*
-            J += (Derivatives.dx(i) * (ElementWeights(i) * dFydy) * Derivatives.dx(i).transpose())
-                 + (Derivatives.dy(i) * (ElementWeights(i) * dFxdx) * Derivatives.dy(i).transpose())
-                 + (Derivatives.dx(i) * (ElementWeights(i) * dFxdy) * Derivatives.dy(i).transpose())
-                 + (Derivatives.dy(i) * (ElementWeights(i) * dFxdy) * Derivatives.dx(i).transpose());
-            */
+            J += Derivatives.dx(i) * (ElementWeights(i) * dFydy).matrix().asDiagonal() * Derivatives.dx(i).transpose()
+                 + Derivatives.dy(i) * (ElementWeights(i) * dFxdx).matrix().asDiagonal() * Derivatives.dy(i).transpose()
+                 + Derivatives.dx(i) * (ElementWeights(i) * dFxdy).matrix().asDiagonal() * Derivatives.dy(i).transpose()
+                 + Derivatives.dy(i) * (ElementWeights(i) * dFxdy).matrix().asDiagonal() * Derivatives.dx(i).transpose();
+
         }
     };
 
@@ -142,7 +143,7 @@ public:
         Derivatives = Domain.template derivative<QuadratureOrder>();
 
         for (size_t i = 0; i != TriangleQuadratureRule<QuadratureOrder>::size; ++i) {
-            ElementWeights(i) = ElementWeights(i) * TriangleQuadratureRule<QuadratureOrder>::w[i];
+            ElementWeights(i) *= TriangleQuadratureRule<QuadratureOrder>::w[i];
         }
     };
 

src/Physics/src/PhysicsConstants.h

+#ifndef OERSTED_PHYSICSCONSTANTS_H
+#define OERSTED_PHYSICSCONSTANTS_H
+
+static constexpr double speed_of_light = 299792458.0;
+static constexpr double mu_0 = 4 * M_PI * 1e-7;
+static constexpr double epsilon_0 = 1.0 / (mu_0 * speed_of_light * speed_of_light);
+
+#endif //OERSTED_PHYSICSCONSTANTS_H

src/Physics/src/Region.h

@@ -3,6 +3,8 @@
 
 #include <vector>
 
+#include "MaterialProperties.h"
+
 template<size_t Dimension>
 class Region {
 };
@@ -10,14 +12,21 @@ class Region {
 template<>
 class Region<2> { // TODO: Rename Triangles to Elements?
 public:
-    Region(std::vector<size_t> tris) : Triangles{tris} {};
+    Region(std::vector<size_t> tris) : Triangles{tris}, Material{Air()} {};
+    Region(std::vector<size_t> tris, MaterialProperties mat) : Triangles{tris}, Material{mat} {};
 
     std::vector<size_t> const &triangles() const { return Triangles; };
 
     size_t const &triangle(size_t i) const { return Triangles[i]; };
 
+    MaterialProperties &material() { return Material; }; // non-const because internal state of the material may be altered
+
+    void material(MaterialProperties mat) { Material = mat; };
+
 protected:
     std::vector<size_t> Triangles;
+
+    MaterialProperties Material;
 };
 
 #endif //OERSTED_REGION_H

src/Physics/src/Triangle.h

@@ -108,7 +108,7 @@ void Triangle<1>::determinant(DiagonalMatrixGroup<TriangleQuadratureRule<Q>::siz
         double xy = p0.y() - p2.y();
         double yx = p1.x() - p2.x();
         double yy = p1.y() - p2.y();
-        mat(i).diagonal()(ID) = xx * yy - xy * yx;
+        mat(i)(ID) = xx * yy - xy * yx;
     }
 }
 

test/Physics/test_FiniteElementMesh.cpp

@@ -34,9 +34,9 @@ TEST_F(MasterTriangle, mass_matrix) {
     auto mat = femesh.basis<2>();
     auto det = femesh.determinant<2>();
 
-    Eigen::SparseMatrix<double> M = (mat[0] * det(0) * mat[0].transpose()) * TriangleQuadratureRule<2>::w[0]
-                                    + (mat[1] * det(1) * mat[1].transpose()) * TriangleQuadratureRule<2>::w[1]
-                                    + (mat[2] * det(2) * mat[2].transpose()) * TriangleQuadratureRule<2>::w[2];
+    Eigen::SparseMatrix<double> M = (mat[0] * det(0).matrix().asDiagonal() * mat[0].transpose()) * TriangleQuadratureRule<2>::w[0]
+                                    + (mat[1] * det(1).matrix().asDiagonal() * mat[1].transpose()) * TriangleQuadratureRule<2>::w[1]
+                                    + (mat[2] * det(2).matrix().asDiagonal() * mat[2].transpose()) * TriangleQuadratureRule<2>::w[2];
 
     EXPECT_NEAR(M.coeffRef(0, 0), 1.0 / 12.0, FLT_EPSILON);
     EXPECT_NEAR(M.coeffRef(1, 0), 1.0 / 24.0, FLT_EPSILON);
@@ -53,7 +53,9 @@ TEST_F(MasterTriangle, stiffness_matrix) {
     auto df = femesh.derivative<1>();
     auto det = femesh.determinant<1>();
 
-    Eigen::SparseMatrix<double> K = (df.dx(0) * det(0) * df.dx(0).transpose() + df.dy(0) * det(0) * df.dy(0).transpose()) * TriangleQuadratureRule<1>::w[0];
+    Eigen::SparseMatrix<double> K =
+            (df.dx(0) * det(0).matrix().asDiagonal() * df.dx(0).transpose() + df.dy(0) * det(0).matrix().asDiagonal() * df.dy(0).transpose()) *
+            TriangleQuadratureRule<1>::w[0];
 
     EXPECT_NEAR(K.coeffRef(0, 0), 1.0, FLT_EPSILON);
     EXPECT_NEAR(K.coeffRef(1, 0), -1.0 / 2.0, FLT_EPSILON);
@@ -83,11 +85,11 @@ TEST_F(MasterTriangle, magnetostatic) {
     auto v = ms.init_unknown_vector();
     auto r = ms.init_unknown_vector();
     auto f = ms.init_unknown_vector();
-    auto Fx = ms.init_element_vector();
-    auto Fy = ms.init_element_vector();
-    auto dFxdx = ms.init_element_matrix();
-    auto dFydy = ms.init_element_matrix();
-    auto dFxdy = ms.init_element_matrix();
+    auto Fx = ms.init_element_array();
+    auto Fy = ms.init_element_array();
+    auto dFxdx = ms.init_element_array();
+    auto dFydy = ms.init_element_array();
+    auto dFxdy = ms.init_element_array();
 
     ms.calculate_forcing(f);
     ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
@@ -507,9 +509,55 @@ public:
     std::vector<Boundary<2>> bnd;
 
     FiniteElementMesh<2, 1> femesh;
+
+    MaterialProperties Linear1000 = MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1000.0));
+
+    void test_flux_density_element_order_1(Eigen::ArrayXd v, Eigen::ArrayXd Bx, Eigen::ArrayXd By, Eigen::ArrayXd Bmag, Eigen::ArrayXd Bang) {
+        double x_q = (1.0 + cos(M_PI / 3.0)) / 3.0;
+        double y_q = sin(M_PI / 3.0) / 3.0;
+        double r_q = hypot(x_q, y_q);
+        double B_q = mu_0 * r_q / 2.0;
+
+        for (size_t i = 0; i != 6; ++i) {
+            EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);
+
+            double ang = 120.0 + 60.0 * i;
+            ang -= 360.0 * (ang > 180.0);
+
+            EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
+        }
+
+        x_q = (1.0 + cos(M_PI / 3.0)) / 2.0;
+        y_q = sin(M_PI / 3.0) / 2.0;
+        r_q = 2.0 - hypot(x_q, y_q);
+        B_q = v(1) / r_q;
+
+        for (size_t i = 6; i < 18; i += 2) {
+            EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);
+
+            double ang = 120.0 + 60.0 * (i - 6) / 2.0;
+            ang -= 360.0 * (ang > 180.0);
+
+            EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
+        }
+
+        x_q = (2.0 * cos(M_PI / 6.0) + 2.0 * cos(M_PI / 6.0 * 3.0)) / 2.0;
+        y_q = (2.0 * sin(M_PI / 6.0) + 2.0 * sin(M_PI / 6.0 * 3.0)) / 2.0;
+        r_q = hypot(x_q, y_q) - 1.0;
+        B_q = v(1) / r_q;
+
+        for (size_t i = 7; i < 18; i += 2) {
+            EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);
+
+            double ang = 120.0 + 60.0 * (i - 6) / 2.0;
+            ang -= 360.0 * (ang > 180.0);
+
+            EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
+        }
+    }
 };
 
-TEST_F(TwoRegionHex, magnetostatic_forcing) {
+TEST_F(TwoRegionHex, magnetostatic_forcing_air) {
     Magnetostatic<2, 1, 1, Variable::A> ms{femesh};
     ms.time(0.0);
 
@@ -532,19 +580,16 @@ TEST_F(TwoRegionHex, magnetostatic_forcing) {
     auto J = ms.init_unknown_matrix();
     auto v = ms.init_unknown_vector();
     auto r = ms.init_unknown_vector();
-    auto Fx = ms.init_element_vector();
-    auto Fy = ms.init_element_vector();
-    auto dFxdx = ms.init_element_matrix();
-    auto dFydy = ms.init_element_matrix();
-    auto dFxdy = ms.init_element_matrix();
+    auto Fx = ms.init_element_array();
+    auto Fy = ms.init_element_array();
+    auto dFxdx = ms.init_element_array();
+    auto dFydy = ms.init_element_array();
+    auto dFxdy = ms.init_element_array();
 
     v.setZero();
 
     ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
 
-    //std::cout << J << std::endl;
-    //std::cout << r << std::endl;
-
     // Test with boundary conditions
     ms.apply_conditions();
 
@@ -554,9 +599,9 @@ TEST_F(TwoRegionHex, magnetostatic_forcing) {
     f = ms.init_unknown_vector();
     Fx = ms.init_element_vector();
     Fy = ms.init_element_vector();
-    dFxdx = ms.init_element_matrix();
-    dFydy = ms.init_element_matrix();
-    dFxdy = ms.init_element_matrix();
+    dFxdx = ms.init_element_array();
+    dFydy = ms.init_element_array();
+    dFxdy = ms.init_element_array();
 
     v.setZero();
 
@@ -567,6 +612,7 @@ TEST_F(TwoRegionHex, magnetostatic_forcing) {
     }
 
     ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
+
     Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
     ldlt.compute(J);
     ASSERT_EQ(ldlt.info(), Eigen::Success);
@@ -579,64 +625,92 @@ TEST_F(TwoRegionHex, magnetostatic_forcing) {
     }
 
     // Fundamental solution is 0.25*r^2 + C in forced region (0)
-    EXPECT_NEAR(v(0), v(1) + 0.25, FLT_EPSILON);
+    EXPECT_NEAR(v(0), v(1) + 0.25 * mu_0, FLT_EPSILON);
     for (size_t i = 1; i != 7; ++i) {
         EXPECT_NEAR(v(1), v(i), FLT_EPSILON);
     }
 
     // Test flux-density values
-    // Values are exact in forced region (0)
-    // Values are approximated in unforced region (1)
-    Eigen::VectorXd Bx = ms.derivatives().dy(0).transpose() * v;
-    Eigen::VectorXd By = -ms.derivatives().dx(0).transpose() * v;
+    Eigen::ArrayXd Bx = ms.derivatives().dy(0).transpose() * v;
+    Eigen::ArrayXd By = -ms.derivatives().dx(0).transpose() * v;
 
-    Eigen::VectorXd Bmag(By.size());
-    Eigen::VectorXd Bang(By.size());
+    Eigen::ArrayXd Bmag(By.size());
+    Eigen::ArrayXd Bang(By.size());
 
     for (size_t i = 0; i != By.size(); ++i) {
         Bmag(i) = hypot(Bx(i), By(i));
         Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
     }
 
-    double x_q = (1.0 + cos(M_PI / 3.0)) / 3.0;
-    double y_q = sin(M_PI / 3.0) / 3.0;
-    double r_q = hypot(x_q, y_q);
-    double B_q = r_q / 2.0;
+    test_flux_density_element_order_1(v, Bx, By, Bmag, Bang);
+}
 
-    for (size_t i = 0; i != 6; ++i) {
-        EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);
+TEST_F(TwoRegionHex, magnetostatic_forcing_1000) {
+    femesh.region(1).material(Linear1000);
 
-        double ang = 120.0 + 60.0 * i;
-        ang -= 360.0 * (ang > 180.0);
+    Magnetostatic<2, 1, 1, Variable::A> ms{femesh};
+    ms.time(0.0);
+    ms.build_matrices();
 
-        EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
-    }
+    auto ff = [](double t) { return 1.0; };
+
+    ms.add_current_density(ff, std::vector<size_t>{0});
+    ms.add_magnetic_insulation(std::vector<size_t>{1});
+
+    ms.apply_conditions();
+
+    auto J = ms.init_unknown_matrix();
+    auto v = ms.init_unknown_vector();
+    auto r = ms.init_unknown_vector();
+    auto f = ms.init_unknown_vector();
+    auto Fx = ms.init_element_array();
+    auto Fy = ms.init_element_array();
+    auto dFxdx = ms.init_element_array();
+    auto dFydy = ms.init_element_array();
+    auto dFxdy = ms.init_element_array();
+
+    v.setZero();
 
-    x_q = (1.0 + cos(M_PI / 3.0)) / 2.0;
-    y_q = sin(M_PI / 3.0) / 2.0;
-    r_q = 2.0 - hypot(x_q, y_q);
-    B_q = v(1) / r_q;
+    ms.calculate_forcing(f);
 
-    for (size_t i = 6; i < 18; i += 2) {
-        EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);
+    ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
 
-        double ang = 120.0 + 60.0 * (i - 6) / 2.0;
-        ang -= 360.0 * (ang > 180.0);
+    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
+    ldlt.compute(J);
+    ASSERT_EQ(ldlt.info(), Eigen::Success);
+    v -= ldlt.solve(r);
 
-        EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
+    // Verify equation is solved
+    ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
+    for (size_t i = 0; i != r.size(); ++i) {
+        EXPECT_NEAR(r(i), 0.0, FLT_EPSILON);
     }
 
-    x_q = (2.0 * cos(M_PI / 6.0) + 2.0 * cos(M_PI / 6.0 * 3.0)) / 2.0;
-    y_q = (2.0 * sin(M_PI / 6.0) + 2.0 * sin(M_PI / 6.0 * 3.0)) / 2.0;
-    r_q = hypot(x_q, y_q) - 1.0;
-    B_q = v(1) / r_q;
+    // Fundamental solution is 0.25*r^2 + C in forced region (0)
+    EXPECT_NEAR(v(0), v(1) + 0.25 * mu_0, FLT_EPSILON);
+    for (size_t i = 1; i != 7; ++i) {
+        EXPECT_NEAR(v(1), v(i), FLT_EPSILON);
+    }
 
-    for (size_t i = 7; i < 18; i += 2) {
-        EXPECT_NEAR(Bmag(i), B_q, FLT_EPSILON);
+    // Test flux-density values
+    Eigen::ArrayXd Bx = ms.derivatives().dy(0).transpose() * v;
+    Eigen::ArrayXd By = -ms.derivatives().dx(0).transpose() * v;
 
-        double ang = 120.0 + 60.0 * (i - 6) / 2.0;
-        ang -= 360.0 * (ang > 180.0);
+    Eigen::ArrayXd Bmag(By.size());
+    Eigen::ArrayXd Bang(By.size());
 
-        EXPECT_NEAR(Bang(i), ang, FLT_EPSILON);
+    for (size_t i = 0; i != By.size(); ++i) {
+        Bmag(i) = hypot(Bx(i), By(i));
+        Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
     }
+
+    // Flux-density in ur=1000 region should be approximately 1000 times greater
+    for (size_t i = 0; i != 6; ++i) {
+        for (size_t j = 6; j != 18; ++j) {
+            EXPECT_GE(Bmag(j), Bmag(i) * 500);
+            EXPECT_LE(Bmag(j), Bmag(i) * 2000);
+        }
+    }
+
+    test_flux_density_element_order_1(v, Bx, By, Bmag, Bang);
 }
\ No newline at end of file