Finished implementation of direct conversion of boundaries from Mesh to FiniteElementMesh

Still need to work on special boundary types (e.g. sliding interfacing and periodic boundaries)

Expanded on magnetostatic testing with a uniform current density on a circular domain

src/Mesh/src/BoundaryConstraint.cpp

@@ -84,4 +84,21 @@ void BoundaryConstraint::refine(Mesh &m, double_t tol) {
 
         m.insert_from_queue();
     }
+}
+
+void BoundaryConstraint::sort_dart_constraints(Mesh const &m) {
+    std::sort(DartConstraints.begin(), DartConstraints.end(),
+              [&](size_t i, size_t j) { return m.dart_constraint(i).S0 < m.dart_constraint(j).S0; }
+    );
+}
+
+std::vector<size_t> BoundaryConstraint::nodes(Mesh const &m) {
+    sort_dart_constraints(m);
+    std::vector<size_t> n;
+    n.reserve(DartConstraints.size());
+    for (size_t i : DartConstraints) {
+        n.push_back(m.dart_constraint(i).base(m));
+    }
+
+    return n;
 }
\ No newline at end of file

src/Mesh/src/BoundaryConstraint.h

@@ -38,8 +38,12 @@ public:
 
     void make_uniform(Mesh &m, double delta_min);
 
+    void sort_dart_constraints(Mesh const &m);
+
     std::shared_ptr<Curve const> curve() const { return ConstraintCurve; };
 
+    std::vector<size_t> nodes(Mesh const &m);
+
 protected:
     std::vector<size_t> DartConstraints;
     std::shared_ptr<Curve const> ConstraintCurve;

src/Mesh/src/DartConstraint.cpp

-#include "DartConstraint.h"
+#include "DartConstraint.h"
\ No newline at end of file
+#include "Mesh.h"
+
+size_t DartConstraint::base(Mesh const &m) const {
+    return (ForwardDart < SIZE_MAX ? m.edge(ForwardDart).node() : m.edge(m.next(ReverseDart)).node());
+};
+
+size_t DartConstraint::tip(Mesh const &m) const {
+    return (ReverseDart < SIZE_MAX ? m.edge(ReverseDart).node() : m.edge(m.next(ForwardDart)).node());
+}
\ No newline at end of file

src/Mesh/src/DartConstraint.h

@@ -18,12 +18,16 @@ public:
 
     size_t self() const { return Self; };
 
-    size_t dart() const { return std::min(ForwardDart, ReverseDart); };
+    size_t dart() const { return (ForwardDart < SIZE_MAX ? ForwardDart : ReverseDart); };
 
     size_t forward_dart() const { return ForwardDart; };
 
     size_t reverse_dart() const { return ReverseDart; };
 
+    size_t base(Mesh const &m) const;
+
+    size_t tip(Mesh const &m) const;
+
     double delta() const { return std::abs(S1 - S0); };
 
     void add_to_queue(Mesh &m) const { Constraint->add_to_queue(m, dart()); };

src/Mesh/src/Edge.cpp

@@ -6,5 +6,5 @@ double_t Edge::length(Mesh const &m) const {
 }
 
 void Edge::add_to_queue(Mesh &m) {
-    m.constraint(Self).add_to_queue(m);
+    m.dart_constraint_from_edge(Self).add_to_queue(m);
 }
\ No newline at end of file

src/Mesh/src/Mesh.cpp

@@ -20,7 +20,7 @@ Mesh::Mesh(Sketch &sketch) {
 
 bool Mesh::are_intersecting(size_t ei, size_t ej) const {
     // TODO, Make more detailed return type enumeration
-    if (is_constrained(ei) && constraint_curve(ei) == constraint_curve(ej)) {
+    if (is_constrained(ei) && curve_from_edge(ei) == curve_from_edge(ej)) {
         return false;
     }
 
@@ -107,7 +107,7 @@ bool Mesh::edges_are_valid() const {
                 result = false;
                 break;
             }
-            if (constraint_curve(e) != constraint_curve(twin(e))) {
+            if (curve_from_edge(e) != curve_from_edge(twin(e))) {
                 result = false;
                 break;
             }

src/Mesh/src/Mesh.h

@@ -86,6 +86,8 @@ public:
 
     size_t size_contours() const { return Contours.size(); };
 
+    size_t size_boundary_constraints() const { return BoundaryConstraints.size(); };
+
     size_t twin(size_t ei) const { return Edges[ei].Twin; };
 
     void add_to_queue(std::unique_ptr<InsertionQueuer> ic) { Queue.push_back(std::move(ic)); };
@@ -95,12 +97,16 @@ public:
     void save_as(std::string path, std::string file_name) const;
 
     // TODO: rename these methods something like, dart_constraint_from_edge_index, curve_from_edge_index
-    DartConstraint const constraint(size_t ei) const { return DartConstraints[Edges[ei].Constraint]; };
+    DartConstraint const dart_constraint_from_edge(size_t ei) const { return DartConstraints[Edges[ei].Constraint]; };
+
+    std::shared_ptr<Curve const> curve_from_edge(size_t ei) const { return DartConstraints[Edges[ei].Constraint].curve(); };
 
-    std::shared_ptr<Curve const> constraint_curve(size_t ei) const { return DartConstraints[Edges[ei].Constraint].curve(); };
+    std::shared_ptr<Curve const> curve(size_t i) const {return BoundaryConstraints[i]->curve(); };
 
     std::shared_ptr<BoundaryConstraint> boundary_constraint(std::shared_ptr<Curve const> const &c) const;
 
+    std::shared_ptr<BoundaryConstraint> boundary_constraint(size_t i) const {return BoundaryConstraints[i]; };
+
     DartConstraint const dart_constraint(size_t i) const { return DartConstraints[i]; };
 
     DartTriangle const &triangle(size_t i) const { return Triangles[i]; };
@@ -138,9 +144,9 @@ public:
 
     Edge const &edge_from_triangle_index(size_t i) const { return Edges[Triangles[i].Edge]; }; // TODO: Rename to triangle_dart
 
-    Edge const &edge_from_triangle(DartTriangle const &dt) { return Edges[dt.Edge]; };
+    Edge const &edge_from_triangle(DartTriangle const &dt) const { return Edges[dt.Edge]; };
 
-    std::array<size_t,3> nodes_from_triangle(DartTriangle const & dt) {
+    std::array<size_t,3> nodes_from_triangle(DartTriangle const & dt) const {
         Edge const &e = edge_from_triangle(dt);
 
         std::array<size_t,3> n;
@@ -152,6 +158,8 @@ public:
         return n;
     };
 
+    std::vector<size_t> boundary_nodes(size_t i) const { return BoundaryConstraints[i]->nodes(*this); }
+
     LocateTriangleResult locate_triangle(Point const p, size_t &ei) const;
 
     LocateTriangleResult locate_triangle(Point const p) const {

src/Physics/src/Boundary.h

@@ -19,14 +19,13 @@ public:
 
     size_t size() const { return Nodes.size(); };
 
-    std::vector<size_t> &nodes() { return Nodes; };
+    auto &nodes() { return Nodes; };
 
-    std::vector<size_t> const &nodes() const { return Nodes; };
+    auto const &nodes() const { return Nodes; };
 
-    size_t node(size_t i) { return Nodes[i]; };
-
-    size_t const &node(size_t i) const { return Nodes[i]; };
+    auto &node(size_t i) { return Nodes[i]; };
 
+    auto const &node(size_t i) const { return Nodes[i]; };
 
 protected:
     std::shared_ptr<Curve const> CurvePtr;

src/Physics/src/FiniteElementMesh.cpp

@@ -27,12 +27,10 @@ template<> FiniteElementMesh<2,1>::FiniteElementMesh(Mesh &m) {
             Then selection can be implemented by passing the shared_ptr that the user has (similar to how the Mesh operations work)
      */
 
-    /*
-    Boundaries.reserve(m.size_constraints());
-    for (size_t i = 0; i!= m.size_constraints(); ++i) {
-        Boundaries.push_back(std::make_shared<Boundary<2>>());
+    Boundaries.reserve(m.size_boundary_constraints());
+    for (size_t i = 0; i!= m.size_boundary_constraints(); ++i) {
+        Boundaries.push_back(std::make_shared<Boundary<2>>(m.curve(i), m.boundary_nodes(i)));
     }
-    */
 
     Regions.reserve(m.size_contours());
     for (size_t i = 0; i != m.size_contours(); ++i) {

src/Physics/src/FiniteElementMesh.h

@@ -26,31 +26,29 @@ public:
 
     size_t size_nodes() const { return Nodes.size(); };
 
-    std::vector<XY> const &nodes() const { return Nodes; };
-
-    XY const &node(size_t i) const { return Nodes[i]; };
-
     size_t size_elements() const { return Triangles.size(); };
 
     std::vector<Triangle<Order>> const &triangles() const { return Triangles; };
 
     Triangle<Order> const &triangle(size_t i) const { return Triangles[i]; };
 
-    auto &regions() { return Regions; };
+    auto &boundaries() { return Boundaries; };
 
-    auto &region(size_t i) { return Regions[i]; };
+    auto &boundary(size_t i) { return Boundaries[i]; };
 
-    auto const &regions() const { return Regions; };
+    auto const &boundaries() const { return Boundaries; };
 
-    auto const &region(size_t i) const { return Regions[i]; };
+    auto const &nodes() const { return Nodes; };
 
-    auto &boundaries() {return Boundaries;};
+    auto const &node(size_t i) const { return Nodes[i]; };
 
-    auto &boundary(size_t i) {return Boundaries[i];};
+    auto &regions() { return Regions; };
 
-    auto const &boundaries() const { return Boundaries; };
+    auto &region(size_t i) { return Regions[i]; };
 
-    auto const &boundary(size_t i) const { return Boundaries[i]; };
+    auto const &regions() const { return Regions; };
+
+    auto const &region(size_t i) const { return Regions[i]; };
 
     template<size_t Q>
     DiagonalMatrixGroup<TriangleQuadrature<Q>::size> determinant() const {
@@ -85,6 +83,18 @@ public:
         return df;
     }
 
+    template<size_t Q>
+    std::vector<std::vector<XY>> quadrature_points() const {
+        std::vector<std::vector<XY>> qp;
+        qp.reserve(Triangles.size());
+
+        for(size_t i = 0; i!= Triangles.size();++i){
+            qp.push_back(Triangles[i].quadrature_points<Q>(Nodes));
+        }
+
+        return qp;
+    }
+
 protected:
     std::vector<XY> Nodes;
     std::vector<Triangle<Order>> Triangles;

src/Physics/src/Physics.h

@@ -213,7 +213,7 @@ public:
                 ++iter;
             }
 
-            // Erase all rows associated with the
+            // Erase all rows associated with the current index
             while (iter != triplets.end() && iter->row() == i) {
                 iter = triplets.erase(iter);
             }
@@ -221,6 +221,11 @@ public:
             ++num_erased;
         }
 
+        while (iter != triplets.end()) { // Finish renumbering the remainder of the entries
+            *iter = Eigen::Triplet<double>(iter->row() - num_erased, iter->col(), iter->value());
+            ++iter;
+        }
+
         size_t cols{Domain.nodes().size()};
         size_t rows{cols - num_erased};
         // TODO: END

src/Physics/src/Triangle.h

@@ -58,6 +58,9 @@ public:
     template<size_t Q>
     void derivative(DerivativeMatrixGroup<TriangleQuadrature<Q>::size> &df, std::vector<XY> const &nodes) const;
 
+    template<size_t Q>
+    std::vector<XY> quadrature_points(std::vector<XY> const& nodes) const;
+
 protected:
     size_t Node[NumNodes];
 };
@@ -144,4 +147,23 @@ void Triangle<1>::derivative(DerivativeMatrixGroup<TriangleQuadrature<Q>::size>
     }
 };
 
+template<>
+template<size_t Q>
+std::vector<XY> Triangle<1>::quadrature_points(std::vector<XY> const& nodes) const {
+    std::vector<XY> qp;
+    qp.reserve(TriangleQuadrature<Q>::size);
+
+    XY const &p0 = nodes[Node[0]];
+    XY const &p1 = nodes[Node[1]];
+    XY const &p2 = nodes[Node[2]];
+
+    for (size_t i = 0; i != TriangleQuadrature<Q>::size; ++i) {
+        double_t x = p0.x() * TriangleQuadrature<Q>::a[i] + p1.x() * TriangleQuadrature<Q>::b[i] + p2.x() * (1.0 - TriangleQuadrature<Q>::a[i] - TriangleQuadrature<Q>::b[i]);
+        double_t y = p0.y() * TriangleQuadrature<Q>::a[i] + p1.y() * TriangleQuadrature<Q>::b[i] + p2.y() * (1.0 - TriangleQuadrature<Q>::a[i] - TriangleQuadrature<Q>::b[i]);
+        qp.emplace_back(x,y);
+    }
+
+    return qp;
+}
+
 #endif //OERSTED_TRIANGLE_H
\ No newline at end of file

test/LibraryIntegration/test_Mesh_to_FEM.cpp

 #include "test_Library_Integration.hpp"
 
-class Library_Integration_Circle : public ::testing::Test {
-public:
-    virtual void SetUp() {
-        sk = Sketch();
+TEST(Library_Integration_Circle, Uniform_Current_Density) {
+    std::string save_dir = "./test/output/Integration/";
 
-        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 * std::cos(M_PI * 2.0 / 3.0), 1.0 * std::sin(M_PI * 2.0 / 3.0));
-        auto v3 = sk.new_element<Vertex>(1.0 * std::cos(-M_PI * 2.0 / 3.0), 1.0 * std::sin(-M_PI * 2.0 / 3.0));
+    // Create Sketch
+    Sketch sk;
 
-        auto c0 = sk.new_element<CircularArc>(v1,v2,v0,1.0);
-        auto c1 = sk.new_element<CircularArc>(v2,v3,v0,1.0);
-        auto c2 = sk.new_element<CircularArc>(v3,v1,v0,1.0);
+    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 * std::cos(M_PI * 2.0 / 3.0), 1.0 * std::sin(M_PI * 2.0 / 3.0));
+    auto v3 = sk.new_element<Vertex>(1.0 * std::cos(-M_PI * 2.0 / 3.0), 1.0 * std::sin(-M_PI * 2.0 / 3.0));
 
-        auto f0 = sk.new_element<Fixation>(v0);
-        auto f1 = sk.new_element<Fixation>(v1);
-        auto f2 = sk.new_element<Fixation>(v2);
-        auto f3 = sk.new_element<Fixation>(v3);
+    auto c0 = sk.new_element<CircularArc>(v1,v2,v0,1.0);
+    auto c1 = sk.new_element<CircularArc>(v2,v3,v0,1.0);
+    auto c2 = sk.new_element<CircularArc>(v3,v1,v0,1.0);
 
-        double_t tol = sk.solve();
-        EXPECT_LE(tol, FLT_EPSILON);
+    auto f0 = sk.new_element<Fixation>(v0);
+    auto f1 = sk.new_element<Fixation>(v1);
+    auto f2 = sk.new_element<Fixation>(v2);
+    auto f3 = sk.new_element<Fixation>(v3);
 
-        bool result = sk.build();
-        EXPECT_TRUE(result);
+    double_t tol = sk.solve();
+    EXPECT_LE(tol, FLT_EPSILON);
 
-        EXPECT_EQ(sk.size_contours(),1);
-    }
+    bool result = sk.build();
+    EXPECT_TRUE(result);
 
-    Sketch sk;
+    EXPECT_EQ(sk.size_contours(),1);
 
-    std::string save_dir = "./test/output/Integration/";
-};
-
-TEST_F(Library_Integration_Circle, Uniform_Current_Density) {
     sk.save_as<SaveMethod::Rasterize>(save_dir, std::string("circle_sketch"));
 
+    // Create Refineable Mesh
     Mesh me{sk};
     me.create();
 
@@ -47,8 +42,97 @@ TEST_F(Library_Integration_Circle, Uniform_Current_Density) {
 
     me.save_as(save_dir, std::string("circle_mesh"));
 
+    // Convert to FiniteElementMesh
+    FiniteElementMesh<2, 1> fem{me};
+
+    for (std::shared_ptr<Boundary<2>> const &b : fem.boundaries() ){
+        double_t a0{-2*M_PI};
+        for(size_t i : b->nodes()) {
+            XY const &p = fem.node(i);
+
+            // Test radius of boundary
+            EXPECT_NEAR(1.0, std::hypot(p.x(), p.y()), FLT_EPSILON);
+
+            // Test ordering of boundary nodes
+            double_t a1{std::atan2(p.y(),p.x())};
+            if (a1 < 0) {
+                a1 += 2 * M_PI;
+            }
+            EXPECT_TRUE(a1 > a0);
+            a0 = a1;
+        }
+    }
+
+    // Create magnetostatic physics
+    Magnetostatic<2, 1, 1, FieldVariable::A> msph{fem};
+
+    msph.add_current_density([](double t) { 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)});
+
+    msph.build_matrices();
+
+    msph.apply_conditions();
+
+    // Initialize matrices
+    auto J = msph.init_unknown_matrix();
+    auto v = msph.init_unknown_vector();
+    auto r = msph.init_unknown_vector();
+    auto f = msph.init_unknown_vector();
+    auto Fx = msph.init_element_array();
+    auto Fy = msph.init_element_array();
+    auto dFxdx = msph.init_element_array();
+    auto dFydy = msph.init_element_array();
+    auto dFxdy = msph.init_element_array();
+
+    //
+    // Set time
+    msph.time(0.0);
+    msph.calculate_forcing(f);
 
-    FiniteElementMesh<2, 1> fem{me}; // TODO: Convert refineable mesh to FiniteElementMesh
+    // Linearize
+    v.setZero();
+    msph.linearize(J, v, r, f, Fx, Fy, dFxdx, dFydy, dFxdy);
 
-    std::cout << std::endl;
+    // Factor
+    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
+    ldlt.compute(J);
+    ASSERT_EQ(ldlt.info(), Eigen::Success);
+
+    // Solve
+    v -= ldlt.solve(r);
+
+    // Test
+    msph.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);
+    }
+
+    // Test flux-density values
+    Eigen::ArrayXd Bx = msph.derivatives().dy(0).transpose() * v;
+    Eigen::ArrayXd By = -msph.derivatives().dx(0).transpose() * v;
+
+    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;
+    }
+
+    std::vector<std::vector<XY>> qp = fem.quadrature_points<1>();
+
+    for (size_t i = 0; i != qp.size(); ++i) {
+        double_t r = std::hypot(qp[i][0].x(), qp[i][0].y());
+
+        EXPECT_NEAR(Bmag(i), r, 0.1);
+
+        double_t a = std::atan2(qp[i][0].y(), qp[i][0].x()) * 180 / M_PI + 90.0;
+
+        if (a > 180.0) { a -= 360.0; }
+        if (a - Bang(i) > 180.0)  { a -= 360.0; }
+        if (a - Bang(i) < -180.0) { a += 360.0; }
+        EXPECT_NEAR(Bang(i), a, 360 * 0.1);
+    }
 }
\ No newline at end of file

test/Mesh/util.cpp

@@ -19,32 +19,32 @@ bool edges_are_valid(Mesh &m) {
 
         if ((e.twin() != e.self())) {
             EXPECT_TRUE(e.node() == m.edge(m.edge(e.twin()).next()).node());
-            EXPECT_TRUE(m.constraint_curve(e.self()) == m.constraint_curve(e.twin()));
-            if (m.constraint_curve(e.self()) != nullptr) {
+            EXPECT_TRUE(m.curve_from_edge(e.self()) == m.curve_from_edge(e.twin()));
+            if (m.curve_from_edge(e.self()) != nullptr) {
                 EXPECT_TRUE(e.orientation() != m.edge(e.twin()).orientation());