Partial implementation of Mesh -> FiniteElementMesh conversion

Nodes, Triangles, Regions are complete (except material properties). Still need to work on Boundaries.

src/Mesh/src/DartTriangle.h

@@ -6,11 +6,16 @@
 
 class DartTriangle {
 public:
-    DartTriangle() : Edge{SIZE_MAX}, Contour{SIZE_MAX} {};
+    DartTriangle() : Contour{SIZE_MAX}, Edge{SIZE_MAX} {};
+
     DartTriangle(size_t e, size_t c) : Edge{e}, Contour{c} {};
 
-    size_t Edge;
+    size_t contour() const { return Contour; };
+
+    size_t edge() const { return Edge; };
+
     size_t Contour;
+    size_t Edge;
 };
 
 #endif //OERSTED_DARTTRIANGLE_H

src/Mesh/src/Mesh.cpp

@@ -796,7 +796,7 @@ void Mesh::mark_triangles() {
 void Mesh::refine_once(std::vector<size_t> index, std::vector<double> radii, std::vector<double> quality) {
     for (size_t i = 0; i != Triangles.size(); ++i) {
         size_t j = index[i];
-        if ((triangle(j).Mark) && ((radii[j] > MaximumElementSize) || (radii[j] > MinimumElementSize && quality[j] < MinimumElementQuality))) {
+        if ((edge_from_triangle_index(j).Mark) && ((radii[j] > MaximumElementSize) || (radii[j] > MinimumElementSize && quality[j] < MinimumElementQuality))) {
             add_circumcenter_to_queue(Triangles[j].Edge);
             insert_from_queue();
         }

src/Mesh/src/Mesh.h

@@ -84,6 +84,8 @@ public:
 
     size_t size_triangles() const { return Triangles.size(); };
 
+    size_t size_contours() const { return Contours.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)); };
@@ -99,10 +101,12 @@ public:
 
     std::shared_ptr<BoundaryConstraint> boundary_constraint(std::shared_ptr<Curve const> const &c) const;
 
-    Point circumcenter(size_t ei) const;
-
     DartConstraint const dart_constraint(size_t i) const { return DartConstraints[i]; };
 
+    DartTriangle const &triangle(size_t i) const { return Triangles[i]; };
+
+    Point circumcenter(size_t ei) const;
+
     Point const &base(Edge const &e) const { return Points[e.Node]; };
 
     Point const &base(size_t i) const { return Points[node(i)]; };
@@ -132,7 +136,21 @@ public:
 
     Edge const &twin(Edge const &e) const { return Edges[e.Twin]; };
 
-    Edge const &triangle(size_t i) const { return Edges[Triangles[i].Edge]; }; // TODO: Rename to triangle_dart
+    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]; };
+
+    std::array<size_t,3> nodes_from_triangle(DartTriangle const & dt) {
+        Edge const &e = edge_from_triangle(dt);
+
+        std::array<size_t,3> n;
+
+        n[0] = e.Node;
+        n[1] = next(e).Node;
+        n[2] = prev(e).Node;
+
+        return n;
+    };
 
     LocateTriangleResult locate_triangle(Point const p, size_t &ei) const;
 

src/Physics/src/Boundary.h

@@ -11,7 +11,7 @@ class Boundary {
 template<>
 class Boundary<2> {
 public:
-    Boundary(std::vector<size_t> nodes) : Nodes{nodes} {};
+    Boundary(std::vector<size_t> &&nodes) : Nodes{nodes} {};
 
     size_t size() const { return Nodes.size(); };
 

src/Physics/src/FiniteElementMesh.cpp

-#include "FiniteElementMesh.h"
\ No newline at end of file
+#include "FiniteElementMesh.h"
+
+template<> FiniteElementMesh<2,1>::FiniteElementMesh(Mesh &m) {
+    // std::vector<XY> nodes
+    // std::vector<Triangle<Order>> tris
+    // std::vector<std::shared_ptr<Region<2>>> r
+    // std::vector<std::shared_ptr<Boundary<2>>> b
+
+    Nodes.reserve(m.size_points()); // TODO: This will be larger when Order > 1
+    for (size_t i = 0; i!= m.size_points(); ++i) {
+        Nodes.emplace_back(m.point(i));
+    }
+
+    /* TODO: Changes to rotational motion implementation
+        Initially, perform a straightforward conversion as currently implemented
+        Define a new Boundary type called (called something like SlidingBoundary)
+            SlidingBoundary holds a reference to the two boundaries and ensures the correct mapping
+        In FiniteElementMesh, add a method called make_sliding(...) or something similar
+            This method will implement a procedure which duplicates the nodes on the boundary and renumbers the nodes in the appropriate elements
+        Then, the SlidingInterface boundary condition holds a pointer to the SlidingBoundary instead of two Boundary objects
+
+        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)
+     */
+
+    /*Boundaries.reserve(m.size_constraints());
+    for (size_t i = 0; i!= m.size_constraionts(); ++i) {
+     ...
+     }
+    */
+    Regions.reserve(m.size_contours());
+    for (size_t i = 0; i != m.size_contours(); ++i) {
+        Regions.push_back(std::make_shared<Region<2>>()); // TODO: Assign material properties here
+    }
+
+    Triangles.reserve(m.size_triangles());
+    for (size_t i = 0; i != m.size_triangles(); ++i) {
+        auto &dt = m.triangle(i);
+        Regions[dt.contour()]->push_back(i);
+        Triangles.emplace_back(i, m.nodes_from_triangle(dt));
+    }
+};
\ No newline at end of file

src/Physics/src/FiniteElementMesh.h

@@ -4,6 +4,8 @@
 #include <memory>
 #include <vector>
 
+#include "Mesh.hpp"
+
 #include "Boundary.h"
 #include "MatrixGroup.h"
 #include "Node.h"
@@ -18,9 +20,9 @@ class FiniteElementMesh<2,Order> {
 public:
     FiniteElementMesh() {};
 
-    FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<Order>> tris,
-                      std::vector<std::shared_ptr<Region<2>>> r, std::vector<std::shared_ptr<Boundary<2>>> b)
-            : Nodes(nodes), Triangles(tris), Regions(r), Boundaries(b) {};
+    FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<Order>> tris, std::vector<std::shared_ptr<Region<2>>> r, std::vector<std::shared_ptr<Boundary<2>>> b) : Nodes(nodes), Triangles(tris), Regions(r), Boundaries(b) {};
+
+    FiniteElementMesh(Mesh &m);
 
     size_t size_nodes() const { return Nodes.size(); };
 

src/Physics/src/Node.h

@@ -5,6 +5,8 @@
 
 #include <cstddef>
 
+#include "Mesh.hpp"
+
 template<size_t N>
 class SmallVector {
 private:
@@ -24,6 +26,8 @@ class SmallVector<2> {
 public:
     SmallVector() : Value{} {};
 
+    SmallVector(Point const& p) : Value{p.X, p.Y} {};
+
     SmallVector(double x, double y) : Value{x, y} {};
 
     double value(size_t i) const { return Value[i]; };

src/Physics/src/Region.h

@@ -12,6 +12,7 @@ class Region {
 template<>
 class Region<2> { // TODO: Rename Triangles to Elements?
 public:
+    Region() : Material{Air()} {};
     Region(std::vector<size_t> elements) : Elements{elements}, Material{Air()} {};
     Region(std::vector<size_t> elements, MaterialProperties material) : Elements{elements}, Material{material} {};
 
@@ -23,6 +24,8 @@ public:
 
     void material(MaterialProperties mat) { Material = mat; };
 
+    void push_back(size_t element) { Elements.push_back(element); };
+
 protected:
     std::vector<size_t> Elements;
 

src/Physics/src/Triangle.h

@@ -38,10 +38,13 @@ public:
 
     Triangle() : Facet{}, Node{} {};
 
-    size_t const &node(size_t const &i) const { return Node[i]; };
+    Triangle(size_t id, std::array<size_t, NumNodes> const &n) : Facet{id} {
+        for (size_t i = 0; i != NumNodes; ++i) {
+            Node[i] = n[i];
+        }
+    };
 
-    template<typename... Args>
-    Triangle(size_t id, Args &&... args) : Facet{id}, Node{std::forward<Args>(args)...} {};
+    size_t const &node(size_t const &i) const { return Node[i]; };
 
     template<size_t D>
     Eigen::Matrix<double, 2, 2> jacobian(std::vector<XY> const &nodes) const; // TODO: Accept Eigen::Matrix<double,2,2> as input

test/LibraryIntegration/test_Mesh_to_FEM.cpp

@@ -48,5 +48,7 @@ TEST_F(Library_Integration_Circle, Uniform_Current_Density) {
     me.save_as(save_dir, std::string("circle_mesh"));
 
 
-    //FiniteElementMesh<2, 1> fem{me}; // TODO: Convert refineable mesh to FiniteElementMesh
+    FiniteElementMesh<2, 1> fem{me}; // TODO: Convert refineable mesh to FiniteElementMesh
+
+    std::cout << std::endl;
 }
\ No newline at end of file

test/Mesh/test_Mesh.cpp

@@ -68,7 +68,7 @@ TEST(Mesh, create_triangle_domain) {
     EXPECT_TRUE(e1.self() == e2.prev());
 
     // Test triangles
-    const Edge e = m.triangle(0);
+    const Edge e = m.edge_from_triangle_index(0);
     Point cc = m.circumcenter(e.self());
     EXPECT_NEAR(0.0, cc.X, TOL);
     EXPECT_NEAR(sqrt(3.0) / 3.0, cc.Y, TOL);
@@ -157,7 +157,7 @@ TEST(Mesh, create_square_domain) {
 
     // Test triangles
     for (size_t i = 0; i < m.size_triangles(); ++i) {
-        Point cc = m.circumcenter(m.triangle(0).self());
+        Point cc = m.circumcenter(m.edge_from_triangle_index(0).self());
         EXPECT_NEAR(0.5, cc.X, TOL);
         EXPECT_NEAR(0.5, cc.Y, TOL);
     }
@@ -229,8 +229,8 @@ TEST(Mesh, create_narrow_diamond_domain) {
     }
 
     // Test triangle circumcenters
-    Point cc0 = m.circumcenter(m.triangle(0).self());
-    Point cc1 = m.circumcenter(m.triangle(1).self());
+    Point cc0 = m.circumcenter(m.edge_from_triangle_index(0).self());
+    Point cc1 = m.circumcenter(m.edge_from_triangle_index(1).self());
 
     EXPECT_NEAR(0.0, cc0.X, TOL);
     EXPECT_NEAR(0.75, std::abs(cc0.Y), TOL);

test/Physics/test_FiniteElementMesh.cpp

 #include "test_Physics.hpp"
 
-inline constexpr size_t operator "" _zu(unsigned long long int i) {
-    return (size_t) i;
-}
-
 class MasterTriangle : public ::testing::Test {
 public:
     virtual void SetUp() {
@@ -11,7 +7,7 @@ public:
         node.push_back(XY{1.0, 0.0});
         node.push_back(XY{0.0, 1.0});
 
-        tri.push_back(Triangle<1>(0, 1_zu, 2_zu, 0_zu));
+        tri.push_back(Triangle<1>(0, {1, 2, 0}));
 
         reg.push_back(std::make_shared<Region<2>>(std::vector<size_t>{0}));
 
@@ -110,13 +106,13 @@ public:
         node.push_back(XY{-1.0, 0.0});
         node.push_back(XY{0.0, -1.0});
 
-        tri.push_back(Triangle<1>(0, 1_zu, 2_zu, 0_zu)); // Master
-        tri.push_back(Triangle<1>(1, 2_zu, 3_zu, 0_zu)); // 90 degree
-        tri.push_back(Triangle<1>(2, 3_zu, 4_zu, 0_zu)); // 180 degree
-        tri.push_back(Triangle<1>(3, 4_zu, 1_zu, 0_zu)); // 270 degree
-        tri.push_back(Triangle<1>(4, 3_zu, 2_zu, 0_zu)); // y-axis reflection
-        tri.push_back(Triangle<1>(5, 1_zu, 4_zu, 0_zu)); // x-axis reflection
-        tri.push_back(Triangle<1>(6, 4_zu, 3_zu, 0_zu)); // x=y reflection
+        tri.push_back(Triangle<1>(0, {1, 2, 0})); // Master
+        tri.push_back(Triangle<1>(1, {2, 3, 0})); // 90 degree
+        tri.push_back(Triangle<1>(2, {3, 4, 0})); // 180 degree
+        tri.push_back(Triangle<1>(3, {4, 1, 0})); // 270 degree
+        tri.push_back(Triangle<1>(4, {3, 2, 0})); // y-axis reflection
+        tri.push_back(Triangle<1>(5, {1, 4, 0})); // x-axis reflection
+        tri.push_back(Triangle<1>(6, {4, 3, 0})); // x=y reflection
 
         reg.push_back(std::make_shared<Region<2>>(std::vector<size_t>{0, 1, 2, 3})); // Triangles 4,5,6 are duplicates
 
@@ -243,8 +239,8 @@ public:
         node.push_back(XY{0.0, 1.0});
         node.push_back(XY{1.0, 1.0});
 
-        tri.push_back(Triangle<1>(0, 0_zu, 1_zu, 2_zu));
-        tri.push_back(Triangle<1>(1, 1_zu, 3_zu, 2_zu));
+        tri.push_back(Triangle<1>(0, {0, 1, 2}));
+        tri.push_back(Triangle<1>(1, {1, 3, 2}));
 
         reg.push_back(std::make_shared<Region<2>>(std::vector<size_t>{0}));
         reg.push_back(std::make_shared<Region<2>>(std::vector<size_t>{1}));
@@ -468,12 +464,12 @@ public:
             node.push_back(XY{r * std::cos(a), r * std::sin(a)});
         }
 
-        tri.push_back(Triangle<1>(0, 0_zu, 1_zu, 2_zu));
-        tri.push_back(Triangle<1>(1, 0_zu, 2_zu, 3_zu));
-        tri.push_back(Triangle<1>(2, 0_zu, 3_zu, 4_zu));
-        tri.push_back(Triangle<1>(3, 0_zu, 4_zu, 5_zu));
-        tri.push_back(Triangle<1>(4, 0_zu, 5_zu, 6_zu));
-        tri.push_back(Triangle<1>(5, 0_zu, 6_zu, 1_zu));
+        tri.push_back(Triangle<1>(0, {0, 1, 2}));
+        tri.push_back(Triangle<1>(1, {0, 2, 3}));
+        tri.push_back(Triangle<1>(2, {0, 3, 4}));
+        tri.push_back(Triangle<1>(3, {0, 4, 5}));
+        tri.push_back(Triangle<1>(4, {0, 5, 6}));
+        tri.push_back(Triangle<1>(5, {0, 6, 1}));
 
         inner_region = std::make_shared<Region<2>>(std::vector<size_t>{0, 1, 2, 3, 4, 5});
         inner_boundary = std::make_shared<Boundary<2>>(std::vector<size_t>{1, 2, 3, 4, 5, 6});
@@ -481,23 +477,23 @@ public:
         reg.push_back(inner_region);
         bnd.push_back(inner_boundary);
 
-        tri.push_back(Triangle<1>(6, 1_zu, 7_zu, 2_zu));
-        tri.push_back(Triangle<1>(7, 2_zu, 7_zu, 8_zu));
+        tri.push_back(Triangle<1>(6, {1, 7, 2}));
+        tri.push_back(Triangle<1>(7, {2, 7, 8}));
 
-        tri.push_back(Triangle<1>(8, 2_zu, 8_zu, 3_zu));
-        tri.push_back(Triangle<1>(9, 3_zu, 8_zu, 9_zu));
+        tri.push_back(Triangle<1>(8, {2, 8, 3}));
+        tri.push_back(Triangle<1>(9, {3, 8, 9}));
 
-        tri.push_back(Triangle<1>(10, 3_zu, 9_zu, 4_zu));
-        tri.push_back(Triangle<1>(11, 4_zu, 9_zu, 10_zu));
+        tri.push_back(Triangle<1>(10, {3, 9, 4}));
+        tri.push_back(Triangle<1>(11, {4, 9, 10}));
 
-        tri.push_back(Triangle<1>(12, 4_zu, 10_zu, 5_zu));
-        tri.push_back(Triangle<1>(13, 5_zu, 10_zu, 11_zu));
+        tri.push_back(Triangle<1>(12, {4, 10, 5}));
+        tri.push_back(Triangle<1>(13, {5, 10, 11}));
 
-        tri.push_back(Triangle<1>(14, 5_zu, 11_zu, 6_zu));
-        tri.push_back(Triangle<1>(15, 6_zu, 11_zu, 12_zu));
+        tri.push_back(Triangle<1>(14, {5, 11, 6}));
+        tri.push_back(Triangle<1>(15, {6, 11, 12}));
 
-        tri.push_back(Triangle<1>(16, 6_zu, 12_zu, 1_zu));
-        tri.push_back(Triangle<1>(17, 1_zu, 12_zu, 7_zu));
+        tri.push_back(Triangle<1>(16, {6, 12, 1}));
+        tri.push_back(Triangle<1>(17, {1, 12, 7}));
 
         outer_region = std::make_shared<Region<2>>(std::vector<size_t>{6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17});
         outer_boundary = std::make_shared<Boundary<2>>(std::vector<size_t>{7, 8, 9, 10, 11, 12});
@@ -741,8 +737,8 @@ public:
         inner_boundary = std::make_shared<Boundary<2>>(std::vector<size_t>{1, 2, 3});
         bnd.push_back(inner_boundary);
 
-        tri.push_back(Triangle<1>(0, 0_zu, 1_zu, 2_zu));
-        tri.push_back(Triangle<1>(1, 0_zu, 2_zu, 3_zu));
+        tri.push_back(Triangle<1>(0, {0, 1, 2}));
+        tri.push_back(Triangle<1>(1, {0, 2, 3}));
 
         inner_region = std::make_shared<Region<2>>(std::vector<size_t>{0, 1});
         reg.push_back(inner_region);
@@ -757,10 +753,10 @@ public:
         outer_boundary = std::make_shared<Boundary<2>>(std::vector<size_t>{4, 5, 6});
         bnd.push_back(outer_boundary);
 
-        tri.push_back(Triangle<1>(2, 1_zu, 4_zu, 2_zu));
-        tri.push_back(Triangle<1>(3, 4_zu, 5_zu, 2_zu));
-        tri.push_back(Triangle<1>(4, 2_zu, 5_zu, 6_zu));
-        tri.push_back(Triangle<1>(5, 2_zu, 6_zu, 3_zu));
+        tri.push_back(Triangle<1>(2, {1, 4, 2}));
+        tri.push_back(Triangle<1>(3, {4, 5, 2}));
+        tri.push_back(Triangle<1>(4, {2, 5, 6}));
+        tri.push_back(Triangle<1>(5, {2, 6, 3}));
 
         outer_region = std::make_shared<Region<2>>(std::vector<size_t>{2, 3, 4, 5});
         reg.push_back(outer_region);
@@ -953,12 +949,12 @@ public:
             node.push_back(XY{r * std::cos(a), r * std::sin(a)});
         }
 
-        tri.push_back(Triangle<1>(0, 0_zu, 1_zu, 2_zu));
-        tri.push_back(Triangle<1>(1, 0_zu, 2_zu, 3_zu));
-        tri.push_back(Triangle<1>(2, 0_zu, 3_zu, 4_zu));
-        tri.push_back(Triangle<1>(3, 0_zu, 4_zu, 5_zu));
-        tri.push_back(Triangle<1>(4, 0_zu, 5_zu, 6_zu));
-        tri.push_back(Triangle<1>(5, 0_zu, 6_zu, 1_zu));
+        tri.push_back(Triangle<1>(0, {0, 1, 2}));
+        tri.push_back(Triangle<1>(1, {0, 2, 3}));
+        tri.push_back(Triangle<1>(2, {0, 3, 4}));
+        tri.push_back(Triangle<1>(3, {0, 4, 5}));
+        tri.push_back(Triangle<1>(4, {0, 5, 6}));
+        tri.push_back(Triangle<1>(5, {0, 6, 1}));
 
         region_ip = std::make_shared<Region<2>>(std::vector<size_t>{0, 1, 2});
         region_in = std::make_shared<Region<2>>(std::vector<size_t>{3, 4, 5});
@@ -981,23 +977,23 @@ public:
             node.push_back(XY{r * std::cos(a), r * std::sin(a)});
         }
 
-        tri.push_back(Triangle<1>(6, 7_zu, 13_zu, 8_zu));
-        tri.push_back(Triangle<1>(7, 8_zu, 13_zu, 14_zu));
+        tri.push_back(Triangle<1>(6, {7, 13, 8}));
+        tri.push_back(Triangle<1>(7, {8, 13, 14}));
 
-        tri.push_back(Triangle<1>(8, 8_zu, 14_zu, 9_zu));
-        tri.push_back(Triangle<1>(9, 9_zu, 14_zu, 15_zu));
+        tri.push_back(Triangle<1>(8, {8, 14, 9}));
+        tri.push_back(Triangle<1>(9, {9, 14, 15}));
 
-        tri.push_back(Triangle<1>(10, 9_zu, 15_zu, 10_zu));
-        tri.push_back(Triangle<1>(11, 10_zu, 15_zu, 16_zu));
+        tri.push_back(Triangle<1>(10, {9, 15, 10}));
+        tri.push_back(Triangle<1>(11, {10, 15, 16}));
 
-        tri.push_back(Triangle<1>(12, 10_zu, 16_zu, 11_zu));
-        tri.push_back(Triangle<1>(13, 11_zu, 16_zu, 17_zu));
+        tri.push_back(Triangle<1>(12, {10, 16, 11}));
+        tri.push_back(Triangle<1>(13, {11, 16, 17}));
 
-        tri.push_back(Triangle<1>(14, 11_zu, 17_zu, 12_zu));
-        tri.push_back(Triangle<1>(15, 12_zu, 17_zu, 18_zu));
+        tri.push_back(Triangle<1>(14, {11, 17, 12}));
+        tri.push_back(Triangle<1>(15, {12, 17, 18}));
 
-        tri.push_back(Triangle<1>(16, 12_zu, 18_zu, 7_zu));
-        tri.push_back(Triangle<1>(17, 7_zu, 18_zu, 13_zu));
+        tri.push_back(Triangle<1>(16, {12, 18, 7}));
+        tri.push_back(Triangle<1>(17, {7, 18, 13}));
 
         region_o = std::make_shared<Region<2>>(std::vector<size_t>{6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17});
         boundary_oi = std::make_shared<Boundary<2>>(std::vector<size_t>{7, 8, 9, 10, 11, 12});