Nightly Commit

Infrastructure and planning for multi-material models

src/Physics/src/FiniteElementMesh.h

@@ -19,6 +19,10 @@ public:
 
     FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<Order>> tris, std::vector<Region<2>> r, std::vector<Boundary<2>> b) : Nodes(nodes), Triangles(tris), Regions(r), Boundaries(b) {};
 
+    size_t size_nodes() const { return Nodes.size(); };
+
+    size_t size_elements() const { return Triangles.size(); };
+
     std::vector<XY> const &nodes() const { return Nodes; };
 
     std::vector<Triangle<Order>> const &triangles() const { return Triangles; };

src/Physics/src/Physics.h

@@ -92,25 +92,47 @@ public:
     using PhysicsData<2, ElementOrder, QuadratureOrder>::weights;
     using PhysicsData<2, ElementOrder, QuadratureOrder>::basis;
 
-    Magnetostatic(FiniteElementMesh<2, ElementOrder> &d) : PhysicsData<2, ElementOrder, QuadratureOrder>{d}, Unknowns{d.nodes().size()} {};
+    Magnetostatic(FiniteElementMesh<2, ElementOrder> &d) : PhysicsData<2, ElementOrder, QuadratureOrder>{d}, Unknowns{d.size_nodes()}, Elements{d.size_elements()} {};
 
-    Eigen::SparseMatrix<double> init_matrix() { return Eigen::SparseMatrix<double>(Unknowns, Unknowns); };
+    Eigen::SparseMatrix<double> init_unknown_matrix() const { return Eigen::SparseMatrix<double>(Unknowns, Unknowns); };
 
-    Eigen::VectorXd init_vector() { return Eigen::VectorXd(Unknowns); };
+    Eigen::DiagonalMatrix<double, Eigen::Dynamic> init_element_matrix() const { return Eigen::DiagonalMatrix<double, Eigen::Dynamic>(Elements); }
 
-    void linearize(Eigen::SparseMatrix<double> &J, Eigen::VectorXd &v, Eigen::VectorXd &r, Eigen::VectorXd &f) {
-        // Calculate residual
-        r = -f;
-        for (size_t i = 0; i!= TriangleQuadratureRule<QuadratureOrder>::size; ++i) {
-            r += Derivatives.dx(i) * (ElementWeights(i) * (Derivatives.dx(i).transpose() * v))
-                 + Derivatives.dy(i) * (ElementWeights(i) * (Derivatives.dy(i).transpose() * v));
-        }
+    Eigen::VectorXd init_unknown_vector() const { return Eigen::VectorXd(Unknowns); };
+
+    Eigen::VectorXd init_element_vector() const { return Eigen::VectorXd(Elements); };
 
-        // Calculate jacobian
+    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) {
+
+        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);
+
+            // Calculate polarization
+            /*
+            for(auto const &r : domain.regions()) {
+                r.material().MagneticFieldIntensity(r.elements(),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 -= ...
+
+            // 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());
+            */
         }
     };
 
@@ -188,6 +210,7 @@ protected:
 
 private:
     size_t Unknowns;
+    size_t Elements;
 };
 
 #endif //OERSTED_PHYSICS_H

src/Physics/src/Region.h

@@ -8,7 +8,7 @@ class Region {
 };
 
 template<>
-class Region<2> {
+class Region<2> { // TODO: Rename Triangles to Elements?
 public:
     Region(std::vector<size_t> tris) : Triangles{tris} {};
 

test/Physics/test_FiniteElementMesh.cpp

@@ -79,13 +79,18 @@ TEST_F(MasterTriangle, magnetostatic) {
 
     ms.build_matrices();
 
-    auto J = ms.init_matrix();
-    auto v = ms.init_vector();
-    auto r = ms.init_vector();
-    auto f = ms.init_vector();
+    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_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();
 
     ms.calculate_forcing(f);
-    ms.linearize(J, v, r, f);
+    ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
 
     //auto ff = [](double t) { return 1.0 * t; };
 
@@ -360,7 +365,7 @@ TEST_F(SimpleSquare, magnetostatic_forcing) {
     ms.build_matrices();
 
     auto ff = [](double t) { return 1.0 * t; };
-    auto f = ms.init_vector();
+    auto f = ms.init_unknown_vector();
 
     ms.add_current_density(ff, std::vector<size_t>{0});
     ms.calculate_forcing(f);
@@ -512,7 +517,7 @@ TEST_F(TwoRegionHex, magnetostatic_forcing) {
     ms.build_matrices();
 
     auto ff = [](double t) { return 1.0; };
-    auto f = ms.init_vector();
+    auto f = ms.init_unknown_vector();
 
     ms.add_current_density(ff, std::vector<size_t>{0});
     ms.add_magnetic_insulation(std::vector<size_t>{1});
@@ -524,12 +529,18 @@ TEST_F(TwoRegionHex, magnetostatic_forcing) {
         EXPECT_NEAR(f(i), nodal_current * 2.0, FLT_EPSILON);
     }
 
-    auto J = ms.init_matrix();
-    auto v = ms.init_vector();
-    auto r = ms.init_vector();
+    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();
+
     v.setZero();
 
-    ms.linearize(J, v, r, f);
+    ms.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
 
     //std::cout << J << std::endl;
     //std::cout << r << std::endl;
@@ -537,10 +548,16 @@ TEST_F(TwoRegionHex, magnetostatic_forcing) {
     // Test with boundary conditions
     ms.apply_conditions();
 
-    J = ms.init_matrix();
-    v = ms.init_vector();
-    r = ms.init_vector();
-    f = ms.init_vector();
+    J = ms.init_unknown_matrix();
+    v = ms.init_unknown_vector();
+    r = ms.init_unknown_vector();
+    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();
+
     v.setZero();
 
     ms.calculate_forcing(f);
@@ -549,14 +566,14 @@ TEST_F(TwoRegionHex, magnetostatic_forcing) {
         EXPECT_NEAR(f(i), nodal_current * 2.0, FLT_EPSILON);
     }
 
-    ms.linearize(J, v, r, f);
+    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);
     v -= ldlt.solve(r);
 
     // Verify equation is solved
-    ms.linearize(J, v, r, f);
+    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);
     }