Store FiniteElementMesh DiscreteBoundary in std::multimap with Curve as key

src/Physics/src/FiniteElementMesh.cpp

@@ -47,9 +47,9 @@ FiniteElementMesh<P, Q>::FiniteElementMesh(Mesh &m) {
     }
     */
 
-    // TODO: CLEAN ME and MESH
+    TODO: CLEAN ME and MESH
 
-    // TODO: LOOK AT WHERE MAP AND SET CAN BE USED INSTEAD OF CRAZY THINGS I HAVE BEEN DOING WITH SHARE_PTR ADDRESS COMPARISONS EG SORT_BOUNDARIES SORT_REGIONS
+    TODO: LOOK AT WHERE MAP AND SET CAN BE USED INSTEAD OF CRAZY THINGS I HAVE BEEN DOING WITH SHARE_PTR ADDRESS COMPARISONS EG SORT_BOUNDARIES SORT_REGIONS
 
     MeshData md = m.mesh_data(P);
 
@@ -58,16 +58,15 @@ FiniteElementMesh<P, Q>::FiniteElementMesh(Mesh &m) {
         Nodes.emplace_back(md.Nodes[i]);
     }
 
-    Boundaries.reserve(md.Boundaries.size());
+    //Boundaries.reserve(md.Boundaries.size());
     for (auto b = ++md.Boundaries.begin(); b != md.Boundaries.end(); ++b) {
-        Boundaries.push_back(std::make_shared<DiscreteBoundary>(b->first, b->second));
+        //Boundaries.push_back(std::make_shared<DiscreteBoundary>(b->first, b->second));
+        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.push_back(std::make_shared<DiscreteRegion>(r->first, r->second)); // TODO: Assign material properties here
         Regions.emplace(r->first, std::make_shared<DiscreteRegion>(r->first, r->second)); // TODO: Assign material properties here
-        // TODO: Emplace
     }
 
     Triangles.resize(md.Edges[0][0].size());
@@ -75,24 +74,12 @@ FiniteElementMesh<P, Q>::FiniteElementMesh(Mesh &m) {
     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].id(k); TODO: FIX THIS STUPID THING
                 Triangles[k].node(n, md.Edges[i][j][k]);
             }
             ++n;
         }
     };
-
-    /*
-    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));
-    }
-    */
-
-    sort_boundaries();
-    //sort_regions();
 };
 
 template<size_t P, size_t Q>
@@ -164,9 +151,11 @@ std::vector<std::shared_ptr<DiscreteBoundary>> FiniteElementMesh<P, Q>::make_dis
 
     // Construct new DiscontinuousBoundary
     bool is_closed = b_vec[0]->is_closed();
-    std::shared_ptr < DiscreteBoundary > db = std::make_shared<DiscreteBoundary>(b_vec[0]->curve(), discontinuous_nodes, is_closed);
-    Boundaries.push_back(db);
-    sort_boundaries();
+    std::shared_ptr<DiscreteBoundary> db = std::make_shared<DiscreteBoundary>(b_vec[0]->curve(), discontinuous_nodes, is_closed);
+    Boundaries.emplace(c, db);
+
+    //Boundaries.push_back(db);
+    //sort_boundaries();
 
     // Renumber nodes in elements when the element is on the "right" side of the plane defined by the discontinuous boundary edge
     // (cross product of edge vector with edge midpoint to incenter vector is negative)
@@ -214,7 +203,9 @@ std::vector<std::shared_ptr<DiscreteBoundary>> FiniteElementMesh<P, Q>::make_dis
 
     // Renumber other boundary nodes when the other boundary is on the right-hand side of the plane defined by the discontinuous boundary edges
     b_vec.push_back(db);
-    for (auto &b : Boundaries) {
+    for (auto iter : Boundaries) {
+        auto b = iter.second;
+
         if (b == b_vec[0] || b == b_vec[1]) {
             continue;
         }

src/Physics/src/FiniteElementMesh.h

@@ -8,7 +8,7 @@ class FiniteElementMesh : public FiniteElementMeshInterface {
 public:
     FiniteElementMesh() {};
 
-    FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<ElementOrder>> tris, std::map<const std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> r, std::vector<std::shared_ptr<DiscreteBoundary>> b) : FiniteElementMeshInterface{nodes,r,b}, Triangles(tris) {};
+    FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<ElementOrder>> tris, std::map<std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> r, std::multimap<std::shared_ptr<const Curve>, std::shared_ptr<DiscreteBoundary>> b) : FiniteElementMeshInterface{nodes,r,b}, Triangles(tris) {};
 
     FiniteElementMesh(Mesh &m);
 

src/Physics/src/FiniteElementMeshInterface.cpp

 #include "FiniteElementMeshInterface.h"
 
-FiniteElementMeshInterface::FiniteElementMeshInterface(std::vector<XY> nodes, std::map<const std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> r, std::vector<std::shared_ptr<DiscreteBoundary>> b)
+FiniteElementMeshInterface::FiniteElementMeshInterface(std::vector<XY> nodes, std::map<std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> r, std::multimap<std::shared_ptr<const Curve>, std::shared_ptr<DiscreteBoundary>> b)
         : Nodes(nodes), Regions(r), Boundaries(b) {
-    sort_boundaries();
-    //sort_regions();
 }
 
 std::vector<std::shared_ptr<DiscreteBoundary>> FiniteElementMeshInterface::boundary(std::shared_ptr<Curve const> const &c) const {
     // Boundary may be discontinuous (e.g. multiple overlapping curves). Therefore, multiple DiscreteBoundaries may be returned (in general, upper != lower)
 
-    auto lower_comp = [](std::shared_ptr<DiscreteBoundary> const &x, size_t const &y) { return (size_t) (x->curve().get()) < y; };
-    auto upper_comp = [](size_t const &y, std::shared_ptr<DiscreteBoundary> const &x) { return y < (size_t) (x->curve().get()); };
+    auto range = Boundaries.equal_range(c);
 
-    auto lower = std::lower_bound(Boundaries.begin(), Boundaries.end(), (size_t) c.get(), lower_comp);
-    auto upper = std::upper_bound(Boundaries.begin(), Boundaries.end(), (size_t) c.get(), upper_comp);
-
-    if (lower != Boundaries.end()) {
-        return std::vector<std::shared_ptr<DiscreteBoundary>>(lower, upper);
-    } else {
-        return std::vector<std::shared_ptr<DiscreteBoundary>>();
-    }
-}
-
-std::shared_ptr<DiscreteRegion> FiniteElementMeshInterface::region(std::shared_ptr<Contour const> const &c) {
-    /*
-    auto comp = [](std::shared_ptr<DiscreteRegion> const &x, size_t const &y) { return (size_t) (x->region().get()) < y; };
-
-    auto iter = std::lower_bound(Regions.begin(), Regions.end(), (size_t) c.get(), comp);
-
-    if (iter != Regions.end()) {
-        return *iter;
-    } else {
-        return std::make_shared<DiscreteRegion>();
+    std::vector<std::shared_ptr<DiscreteBoundary>> output;
+    for (auto it = range.first; it != range.second; ++it) {
+        output.push_back(it->second);
     }
-     */
-    return Regions[c];
-}
 
-void FiniteElementMeshInterface::sort_boundaries() {
-    auto comp = [](std::shared_ptr<DiscreteBoundary> const &x, std::shared_ptr<DiscreteBoundary> const &y) { return (size_t) (x->curve().get()) < (size_t) (y->curve().get()); };
-    std::sort(Boundaries.begin(), Boundaries.end(), comp);
-}
-
-/*
-void FiniteElementMeshInterface::sort_regions() {
-    auto comp = [](std::shared_ptr<DiscreteRegion> const &x, std::shared_ptr<DiscreteRegion> const &y) { return (size_t) (x->region().get()) < (size_t) (y->region().get()); };
-    std::sort(Regions.begin(), Regions.end(), comp);
-}
- */
\ No newline at end of file
+    return output;
+}
\ No newline at end of file

src/Physics/src/FiniteElementMeshInterface.h

@@ -16,13 +16,13 @@ class FiniteElementMeshInterface {
 public:
     FiniteElementMeshInterface() {};
 
-    FiniteElementMeshInterface(std::vector<XY> nodes, std::map<const std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> r, std::vector<std::shared_ptr<DiscreteBoundary>> b);
+    FiniteElementMeshInterface(std::vector<XY> nodes, std::map<std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> r, std::multimap<std::shared_ptr<const Curve>, std::shared_ptr<DiscreteBoundary>> b);
 
     virtual ~FiniteElementMeshInterface() {};
 
     auto &boundaries() { return Boundaries; };
 
-    auto &boundary(size_t i) { return Boundaries[i]; };
+    //auto &boundary(size_t i) { return Boundaries[i]; };
 
     std::vector<std::shared_ptr<DiscreteBoundary>> boundary(std::shared_ptr<Curve const> const &c) const;
 
@@ -40,7 +40,7 @@ public:
 
     auto const &regions() const { return Regions; };
 
-    std::shared_ptr<DiscreteRegion> region(std::shared_ptr<Contour const> const &c);
+    std::shared_ptr<DiscreteRegion> region(std::shared_ptr<Contour const> const &c) { return Regions[c]; };
 
     DiagonalMatrixGroup diagonal_matrix_group() const { return DiagonalMatrixGroup{quadrature_size(), elements_size()}; };
 
@@ -77,17 +77,11 @@ public:
 protected:
     std::vector<XY> Nodes;
 
-    // TODO: USE SET OR MAP TO STORE REGIONS AND BOUNDARIES BASED ON CONTOUR AND CURVE SHARED_PTRS
+    TODO: Define type alias for mapping
 
-    //std::vector<std::shared_ptr<DiscreteRegion>> Regions; // Contains vector of size_t referencing Triangles (and later Quadrilaterals)
+    std::map<std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> Regions;
 
-    std::map<const std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> Regions;
-
-    std::vector<std::shared_ptr<DiscreteBoundary>> Boundaries;
-
-    void sort_boundaries();
-
-    // void sort_regions();
+    std::multimap<std::shared_ptr<const Curve>, std::shared_ptr<DiscreteBoundary>> Boundaries;
 };
 
 #endif //OERSTED_FINITEELEMENTMESHINTERFACE_H

test/LibraryIntegration/Mesh_To_FEM_Test.cpp

@@ -45,7 +45,9 @@ TEST(Full_Circle, Magnetostatic_Uniform_Current_Density) {
     // Convert to FiniteElementMesh
     auto fem = std::make_shared<FiniteElementMesh<1,1>>(me);
 
-    for (std::shared_ptr<DiscreteBoundary> const &b : fem->boundaries()) {
+    for (auto iter : fem->boundaries()) {
+        auto b = iter.second;
+
         double_t a0{-2 * M_PI};
         for (size_t i : b->nodes()) {
             XY const &p = fem->node(i);
@@ -72,7 +74,9 @@ TEST(Full_Circle, Magnetostatic_Uniform_Current_Density) {
     auto cdr = me.select_contour(Point{0.0,0.0});
     msph.add_current_density([](FunctionArguments args) { return 1.0 / (2.0 * M_PI * 1e-7); }, {fem->region(cdr)});
 
-    msph.add_magnetic_insulation({fem->boundary(0), fem->boundary(1), fem->boundary(2)});
+    msph.add_magnetic_insulation(fem->boundary(c0));
+    msph.add_magnetic_insulation(fem->boundary(c1));
+    msph.add_magnetic_insulation(fem->boundary(c2));
 
     msph.build_matrices();
 

test/Physics/Finite_Element_Mesh_Test.cpp

@@ -10,20 +10,24 @@ public:
 
         tri.push_back(Triangle<1>(0, {1, 2, 0}));
 
-        reg.emplace(c[0], std::make_shared<DiscreteRegion>(std::vector<size_t>{0}));
+        reg.emplace(std::make_shared<Contour>(), std::make_shared<DiscreteRegion>(std::vector<size_t>{0}));
 
-        bnd.push_back(std::make_shared<DiscreteBoundary>(std::vector<size_t>{0, 1}));
-        bnd.push_back(std::make_shared<DiscreteBoundary>(std::vector<size_t>{0, 2}));
-        bnd.push_back(std::make_shared<DiscreteBoundary>(std::vector<size_t>{1, 2}));
+        for (auto iter : reg) {
+            creg.push_back(iter.first);
+        }
+
+        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}));
+        bnd.emplace(std::make_shared<LineSegment>(), std::make_shared<DiscreteBoundary>(std::vector<size_t>{1, 2}));
     }
 
     std::vector<XY> node;
     std::vector<Triangle<1>> tri;
 
-    std::vector<std::shared_ptr<const Contour>> c{std::make_shared<const Contour>()};
+    std::vector<std::shared_ptr<const Contour>> creg;
 
-    std::map<const std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> reg;
-    std::vector<std::shared_ptr<DiscreteBoundary>> bnd;
+    std::map<std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> reg;
+    std::multimap<std::shared_ptr<const Curve>, std::shared_ptr<DiscreteBoundary>> bnd;
 
     template<size_t QuadratureOrder>
     std::shared_ptr<FiniteElementMesh<1, QuadratureOrder>> get_mesh() {
@@ -83,7 +87,7 @@ TEST_F(MasterTriangle, magnetostatic) {
 
     auto ff = [](FunctionArguments args) { return 1.0 * args["t"]; };
 
-    ms.add_current_density(ff, {reg[c[0]]});
+    ms.add_current_density(ff, {reg[creg[0]]});
 
     ms.build_matrices();
 
@@ -117,16 +121,14 @@ public:
         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.emplace(c[0], std::make_shared<DiscreteRegion>(std::vector<size_t>{0, 1, 2, 3})); // Triangles 4,5,6 are duplicates
+        reg.emplace(std::make_shared<Contour>(), std::make_shared<DiscreteRegion>(std::vector<size_t>{0, 1, 2, 3})); // Triangles 4,5,6 are duplicates
     }
 
     std::vector<XY> node;
     std::vector<Triangle<1>> tri;
 
-    std::vector<std::shared_ptr<const Contour>> c{std::make_shared<const Contour>()};
-
-    std::map<const std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> reg;
-    std::vector<std::shared_ptr<DiscreteBoundary>> bnd;
+    std::map<std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> reg;
+    std::multimap<std::shared_ptr<const Curve>, std::shared_ptr<DiscreteBoundary>> bnd;
 
     template<size_t QuadratureOrder>
     std::shared_ptr<FiniteElementMesh<1, QuadratureOrder>> get_mesh() {
@@ -248,30 +250,34 @@ public:
         tri.push_back(Triangle<1>(0, {0, 1, 2}));
         tri.push_back(Triangle<1>(1, {1, 3, 2}));
 
-        reg.emplace(c[0], std::make_shared<DiscreteRegion>(std::vector<size_t>{0}));
-        reg.emplace(c[1], std::make_shared<DiscreteRegion>(std::vector<size_t>{1}));
-
-        EXPECT_EQ(reg.size(), 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);
+        }
 
-        bnd.push_back(std::make_shared<DiscreteBoundary>(std::vector<size_t>{0, 1}));
-        bnd.push_back(std::make_shared<DiscreteBoundary>(std::vector<size_t>{0, 2}));
-        bnd.push_back(std::make_shared<DiscreteBoundary>(std::vector<size_t>{1, 2}));
-        bnd.push_back(std::make_shared<DiscreteBoundary>(std::vector<size_t>{1, 3}));
-        bnd.push_back(std::make_shared<DiscreteBoundary>(std::vector<size_t>{2, 3}));
+        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}));
+        bnd.emplace(std::make_shared<LineSegment>(), std::make_shared<DiscreteBoundary>(std::vector<size_t>{1, 2}));
+        bnd.emplace(std::make_shared<LineSegment>(), std::make_shared<DiscreteBoundary>(std::vector<size_t>{1, 3}));
+        bnd.emplace(std::make_shared<LineSegment>(), std::make_shared<DiscreteBoundary>(std::vector<size_t>{2, 3}));
 
         vx_m << 1.0, 0.0, 1.0, 0.0;
         vx_p << 0.0, 1.0, 0.0, 1.0;
         vy_m << 1.0, 1.0, 0.0, 0.0;
         vy_p << 0.0, 0.0, 1.0, 1.0;
+
+        EXPECT_EQ(reg.size(), 2);
+        EXPECT_EQ(bnd.size(), 5);
     }
 
     std::vector<XY> node;
     std::vector<Triangle<1>> tri;
 
-    std::vector<std::shared_ptr<const Contour>> c{std::make_shared<const Contour>(), std::make_shared<const Contour>()};
+    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::map<const std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> reg;
-    std::vector<std::shared_ptr<DiscreteBoundary>> bnd;
+    std::vector<std::shared_ptr<const Contour>> creg;
 
     Eigen::Vector4d vx_m;
     Eigen::Vector4d vx_p;
@@ -380,7 +386,7 @@ TEST_F(SimpleSquare, magnetostatic_forcing) {
     auto ff = [](FunctionArguments args) { return 1.0 * args["t"]; };
     auto f = ms.init_unknown_vector();
 
-    ms.add_current_density(ff, {reg[c[0]]});
+    ms.add_current_density(ff, {reg[creg[0]]});
     ms.calculate_forcing(f, fargs);
     EXPECT_NEAR(f(0), 1.0 / 6.0, FLT_EPSILON);
     EXPECT_NEAR(f(1), 1.0 / 6.0, FLT_EPSILON);
@@ -388,7 +394,7 @@ TEST_F(SimpleSquare, magnetostatic_forcing) {
     EXPECT_NEAR(f(3), 0.0 / 6.0, FLT_EPSILON);
     ms.remove_current_density(0);
 
-    ms.add_current_density(ff, {reg[c[1]]});
+    ms.add_current_density(ff, {reg[creg[1]]});
     ms.calculate_forcing(f, fargs);
     EXPECT_NEAR(f(0), 0.0 / 6.0, FLT_EPSILON);
     EXPECT_NEAR(f(1), 1.0 / 6.0, FLT_EPSILON);
@@ -396,7 +402,7 @@ TEST_F(SimpleSquare, magnetostatic_forcing) {
     EXPECT_NEAR(f(3), 1.0 / 6.0, FLT_EPSILON);
     ms.remove_current_density(0);
 
-    ms.add_current_density(ff, {reg[c[0]], reg[c[1]]});
+    ms.add_current_density(ff, {reg[creg[0]], reg[creg[1]]});
     ms.calculate_forcing(f, fargs);
     EXPECT_NEAR(f(0), 1.0 / 6.0, FLT_EPSILON);
     EXPECT_NEAR(f(1), 2.0 / 6.0, FLT_EPSILON);
@@ -464,6 +470,14 @@ TEST_F(SimpleSquare, magnetostatic_forcing) {
 class TwoRegionHexagon : public ::testing::Test {
 public:
     virtual void SetUp() {
+        for (size_t i = 0; i != 2; ++i) {
+            creg.push_back(std::make_shared<Contour>());
+        }
+
+        for (size_t i = 0; i != 2; ++i) {
+            cbnd.push_back(std::make_shared<LineSegment>());
+        }
+
         node.push_back(XY{0.0, 0.0});
 
         double r{1.0};
@@ -488,8 +502,8 @@ public:
         inner_region = std::make_shared<DiscreteRegion>(std::vector<size_t>{0, 1, 2, 3, 4, 5});
         inner_boundary = std::make_shared<DiscreteBoundary>(std::vector<size_t>{1, 2, 3, 4, 5, 6});
 
-        reg.emplace(c[0], inner_region);
-        bnd.push_back(inner_boundary);
+        reg.emplace(creg[0], inner_region);
+        bnd.emplace(cbnd[0], inner_boundary);
 
         tri.push_back(Triangle<1>(6, {1, 7, 2}));
         tri.push_back(Triangle<1>(7, {2, 7, 8}));
@@ -512,10 +526,11 @@ public:
         outer_region = std::make_shared<DiscreteRegion>(std::vector<size_t>{6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17});
         outer_boundary = std::make_shared<DiscreteBoundary>(std::vector<size_t>{7, 8, 9, 10, 11, 12});
 
-        reg.emplace(c[1], outer_region);
-        bnd.push_back(outer_boundary);
+        reg.emplace(creg[1], outer_region);
+        bnd.emplace(cbnd[1], outer_boundary);
 
         EXPECT_EQ(reg.size(), 2);
+        EXPECT_EQ(bnd.size(), 2);
 
         femesh = std::make_shared<FiniteElementMesh<1, 1>>(node, tri, reg, bnd);
     }
@@ -523,10 +538,11 @@ public:
     std::vector<XY> node;
     std::vector<Triangle<1>> tri;
 
-    std::vector<std::shared_ptr<const Contour>> c{std::make_shared<const Contour>(),std::make_shared<const Contour>()};
+    std::vector<std::shared_ptr<const Contour>> creg;
+    std::vector<std::shared_ptr<const Curve>> cbnd;
 
-    std::map<const std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> reg;
-    std::vector<std::shared_ptr<DiscreteBoundary>> bnd;
+    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::shared_ptr<DiscreteRegion> inner_region;
     std::shared_ptr<DiscreteRegion> outer_region;
@@ -737,6 +753,14 @@ TEST_F(TwoRegionHexagon, magnetostatic_forcing_1000) {
 class Two_Region_Quarter_Octagon_Continuous : public ::testing::Test {
 public:
     virtual void SetUp() {
+        for (size_t i = 0; i != 2; ++i) {
+            creg.push_back(std::make_shared<Contour>());
+        }
+
+        for (size_t i = 0; i != 6; ++i) {
+            cbnd.push_back(std::make_shared<LineSegment>());
+        }
+
         // Inner Region
         node.push_back(XY{0.0, 0.0});
 
@@ -747,13 +771,13 @@ public:
         }
 
         inner_boundary = std::make_shared<DiscreteBoundary>(std::vector<size_t>{1, 2, 3});
-        bnd.push_back(inner_boundary);
+        bnd.emplace(cbnd[0], inner_boundary);
 
         tri.push_back(Triangle<1>(0, {0, 1, 2}));
         tri.push_back(Triangle<1>(1, {0, 2, 3}));
 
         inner_region = std::make_shared<DiscreteRegion>(std::vector<size_t>{0, 1});
-        reg.emplace(c[0], inner_region);
+        reg.emplace(creg[0], inner_region);
 
         // Outer Region
         r = 2.0;
@@ -763,7 +787,7 @@ public:
         }
 
         outer_boundary = std::make_shared<DiscreteBoundary>(std::vector<size_t>{4, 5, 6});
-        bnd.push_back(outer_boundary);
+        bnd.emplace(cbnd[1], outer_boundary);
 
         tri.push_back(Triangle<1>(2, {1, 4, 2}));
         tri.push_back(Triangle<1>(3, {4, 5, 2}));
@@ -771,7 +795,7 @@ public:
         tri.push_back(Triangle<1>(5, {2, 6, 3}));
 
         outer_region = std::make_shared<DiscreteRegion>(std::vector<size_t>{2, 3, 4, 5});
-        reg.emplace(c[1], outer_region);
+        reg.emplace(creg[1], outer_region);
 
         // Periodic Boundary
         pbi0 = std::make_shared<DiscreteBoundary>(std::vector<size_t>{0,1});
@@ -779,23 +803,25 @@ public:
         pbo0 = std::make_shared<DiscreteBoundary>(std::vector<size_t>{1,4});
         pbo1 = std::make_shared<DiscreteBoundary>(std::vector<size_t>{3,6});
 
-        bnd.push_back(pbi0);
-        bnd.push_back(pbi1);
-        bnd.push_back(pbo0);
-        bnd.push_back(pbo1);
+        bnd.emplace(cbnd[2],pbi0);
+        bnd.emplace(cbnd[3],pbi1);
+        bnd.emplace(cbnd[4],pbo0);
+        bnd.emplace(cbnd[5],pbo1);
 
         femesh = std::make_shared<FiniteElementMesh<1, 1>>(node, tri, reg, bnd);
 
         EXPECT_EQ(reg.size(), 2);
+        EXPECT_EQ(bnd.size(), 6);
     }
 
     std::vector<XY> node;
     std::vector<Triangle<1>> tri;
 
-    std::vector<std::shared_ptr<const Contour>> c{std::make_shared<const Contour>(),std::make_shared<const Contour>()};
+    std::vector<std::shared_ptr<const Contour>> creg;
+    std::vector<std::shared_ptr<const Curve>> cbnd;
 
-    std::map<const std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> reg;
-    std::vector<std::shared_ptr<DiscreteBoundary>> bnd;
+    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::shared_ptr<DiscreteRegion> inner_region;
     std::shared_ptr<DiscreteRegion> outer_region;
@@ -960,6 +986,14 @@ TEST_F(Two_Region_Quarter_Octagon_Continuous, magnetostatic_antiperiodic) {
 class Two_Region_Hexagon_Rotation_Discontinuous : public ::testing::Test {
 public:
     virtual void SetUp() {
+        for (size_t i = 0; i != 3; ++i) {
+            creg.push_back(std::make_shared<Contour>());
+        }
+
+        for (size_t i = 0; i != 3; ++i) {
+            cbnd.push_back(std::make_shared<LineSegment>());
+        }
+
         // Inner Region
         node.push_back(XY{0.0, 0.0});
 
@@ -977,13 +1011,13 @@ public:
         tri.push_back(Triangle<1>(5, {0, 6, 1}));
 
         region_ip = std::make_shared<DiscreteRegion>(std::vector<size_t>{0, 1, 2});
-        reg.emplace(c[0], region_ip);
+        reg.emplace(creg[0], region_ip);
 
         region_in = std::make_shared<DiscreteRegion>(std::vector<size_t>{3, 4, 5});
-        reg.emplace(c[1], region_in);
+        reg.emplace(creg[1], region_in);
 
         inner_boundary_0 = std::make_shared<DiscreteBoundary>(std::vector<size_t>{1, 2, 3, 4, 5, 6, 1}, true);
-        bnd.push_back(inner_boundary_0);
+        bnd.emplace(cbnd[0], inner_boundary_0);
 
         // Discontinuous Boundary
         r = 1.0;
@@ -993,7 +1027,7 @@ public:
         }
 
         inner_boundary_1 = std::make_shared<DiscreteBoundary>(std::vector<size_t>{7, 8, 9, 10, 11, 12, 7}, true);
-        bnd.push_back(inner_boundary_1);
+        bnd.emplace(cbnd[1], inner_boundary_1);
 
         // Outer Region
         r = 2.0;
@@ -1021,23 +1055,25 @@ public:
         tri.push_back(Triangle<1>(17, {7, 18, 13}));
 
         region_o = std::make_shared<DiscreteRegion>(std::vector<size_t>{6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17});
-        reg.emplace(c[2], region_o);
+        reg.emplace(creg[2], region_o);
 
         outer_boundary = std::make_shared<DiscreteBoundary>(std::vector<size_t>{13, 14, 15, 16, 17, 18}, true);
-        bnd.push_back(outer_boundary);
+        bnd.emplace(cbnd[2], outer_boundary);
 
         femesh = std::make_shared<FiniteElementMesh<1, 1>>(node, tri, reg, bnd);
 
         EXPECT_EQ(reg.size(), 3);
+        EXPECT_EQ(bnd.size(), 3);
     }
 
     std::vector<XY> node;
     std::vector<Triangle<1>> tri;
 
-    std::vector<std::shared_ptr<const Contour>> c{std::make_shared<const Contour>(),std::make_shared<const Contour>(),std::make_shared<const Contour>()};
+    std::vector<std::shared_ptr<const Contour>> creg;
+    std::vector<std::shared_ptr<const Curve>> cbnd;
 
-    std::map<const std::shared_ptr<const Contour>, std::shared_ptr<DiscreteRegion>> reg;
-    std::vector<std::shared_ptr<DiscreteBoundary>> bnd;
+    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::shared_ptr<DiscreteRegion> region_ip;
     std::shared_ptr<DiscreteRegion> region_in;
@@ -1193,6 +1229,14 @@ TEST_F(Two_Region_Hexagon_Rotation_Discontinuous, magnetostatic_forcing) {
 class Two_Region_Quarter_Octagon_Discontinuous : public ::testing::Test {
 public:
     virtual void SetUp() {
+        for (size_t i = 0; i != 2; ++i) {
+            creg.push_back(std::make_shared<Contour>());
+        }
+
+        for (size_t i = 0; i != 6; ++i) {
+            cbnd.push_back(std::make_shared<LineSegment>());
+        }
+
         // Inner Region
         node.push_back(XY{0.0, 0.0});
 
@@ -1203,13 +1247,13 @@ public:
         }