Preparation for 2nd order elements

Store FiniteElementMesh DiscreteRegions in std::map with Contour as key

src/Mesh/src/Mesh.h

@@ -27,6 +27,38 @@ enum class InsertPointResult {
     Success, Midpoint, Duplicate, Failure
 };
 
+class NodeData {
+public:
+    NodeData(){};
+
+    NodeData(size_t elid, size_t edid, size_t noid, std::shared_ptr<const Curve> cc, double_t s, Point p) : ElementID{elid}, EdgeID{edid}, NodeID{noid}, Constraint{cc}, Parameter{s}, Coordinates{p} {};
+
+    friend bool operator<(const NodeData &en0, const NodeData &en1) { return (en0.Coordinates < en1.Coordinates); };
+
+    size_t ElementID;
+
+    size_t EdgeID;
+
+    size_t NodeID;
+
+    std::shared_ptr<const Curve> Constraint;
+
+    double_t Parameter;
+
+    Point Coordinates;
+};
+
+class MeshData {
+public:
+    std::vector<Point> Nodes;
+
+    std::vector<std::vector<std::vector<size_t>>> Edges;
+
+    std::map<std::shared_ptr<const Curve>, std::set<std::pair<double_t,size_t>>> Boundaries;
+
+    std::map<std::shared_ptr<const Contour>, std::vector<size_t>> Regions;
+};
+
 class Mesh { // TODO: Namespaces. Also, need to rename a bunch of things to be more consistent and clear w.r.t. how various arrays are accessed (e.g. by edge/dart index versus plain index)
 public:
     Mesh() {};
@@ -115,6 +147,8 @@ public:
 
     std::shared_ptr<Contour const> contour(size_t i) const { return Contours[i]; };
 
+    std::shared_ptr<Contour const> contour_from_triangle_index(size_t i) const { return Contours[Triangles[i].Contour]; };
+
     std::shared_ptr<Curve const> curve_from_edge(size_t ei) const { return DartConstraints[Edges[ei].Constraint].curve(); };
 
     std::shared_ptr<Curve const> curve(size_t i) const {return BoundaryConstraints[i]->curve(); };
@@ -199,6 +233,78 @@ public:
         return n;
     };
 
+    MeshData mesh_data(size_t p) const {
+        std::vector<NodeData> data;
+        std::vector<Point> points;
+        std::vector<std::vector<std::vector<size_t>>> edges;
+
+        std::map<std::shared_ptr<const Curve>, std::set<std::pair<double_t,size_t>>> boundaries;
+        std::map<std::shared_ptr<const Contour>, std::vector<size_t>> regions;
+
+        for (size_t t = 0; t != size_triangles(); ++t) {
+            Edge e = edge_from_triangle_index(t);
+            regions[contour_from_triangle_index(t)].push_back(t);
+
+            for (size_t k = 0; k != 3; ++k) {
+                DartConstraint dc = dart_constraint_from_edge(e.Self);
+                Point point0 = base(e);
+                Point point1 = tip(e);
+
+                double_t param0, param1;
+                if (e.orientation()) {
+                    param0 = dc.S0;
+                    param1 = dc.S1;
+                } else {
+                    param0 = dc.S1;
+                    param1 = dc.S0;
+                }
+
+                for (size_t i = 0; i != p + 1; i++) {
+                    double_t s = i / (double_t) (p);
+                    Point pointi = point0 * s + point1 * (1.0 - s);
+                    double_t parami = param0 * s + param1 * (1.0 - s);
+
+                    data.emplace_back(t, k, i, dc.curve(), parami, pointi);
+                }
+
+                e = next(e);
+            }
+        }
+
+        std::sort(data.begin(), data.end());
+
+        edges.resize(p + 1);
+        for (size_t i = 0; i != p + 1; ++i) {
+            edges[i].resize(3);
+            for (size_t j = 0; j != 3; ++j) {
+                edges[i][j].resize(size_triangles());
+            }
+        }
+
+        size_t i = 0;
+        size_t n = 0;
+        while (i < data.size()) {
+            points.push_back(data[i].Coordinates);
+
+            edges[data[i].NodeID][data[i].EdgeID][data[i].ElementID] = n;
+            boundaries[data[i].Constraint].emplace(data[i].Parameter,n);
+
+            size_t j = i + 1;
+            while (j < data.size() && data[i].Coordinates == data[j].Coordinates) {
+                edges[data[j].NodeID][data[j].EdgeID][data[j].ElementID] = n;
+                boundaries[data[j].Constraint].emplace(data[j].Parameter,n);
+
+                ++j;
+            }
+
+            i = j;
+            ++n;
+        }
+        edges.erase(edges.end()-1);
+
+        return MeshData{points,edges,boundaries,regions};
+    }
+
     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;

src/Mesh/src/Point.h

@@ -44,6 +44,13 @@ public:
         return *this;
     }
 
+    Point &operator*=(double rhs) {
+        X *= rhs;
+        Y *= rhs;
+
+        return *this;
+    }
+
     friend Point operator+(Point lhs, Point const &rhs) {
         lhs += rhs;
         return lhs;
@@ -53,6 +60,15 @@ public:
         p /= den;
         return p;
     }
+
+    friend Point operator*(Point p, double const &rhs) {
+        p *= rhs;
+        return p;
+    }
+
+    friend bool operator<(const Point &p0, const Point &p1) {
+        return ((p0.X < p1.X)) || ((p0.X == p1.X) && (p0.Y < p1.Y));
+    }
 };
 
 double dist(Point const &p0, Point const &p1);

src/Physics/src/DiscreteBoundary.h

@@ -2,6 +2,7 @@
 #define OERSTED_BOUNDARY_H
 
 #include <vector>
+#include <set>
 #include <cstddef>
 
 #include "Sketch.hpp"
@@ -14,7 +15,12 @@ public:
 
     DiscreteBoundary(std::shared_ptr<Curve const> cptr, std::vector<size_t> nodes, bool closed = false) : CurvePtr{cptr}, Nodes{nodes}, Closed{closed} {};
 
-    virtual ~DiscreteBoundary() {};
+    DiscreteBoundary(std::shared_ptr<const Curve> cptr, std::set<std::pair<double_t, size_t>> nodes, bool closed = false) : CurvePtr{cptr}, Closed{closed} {
+        Nodes.reserve(nodes.size());
+        for (auto it = nodes.begin(); it != nodes.end(); ++it) {
+            Nodes.push_back(it->second);
+        }
+    }
 
     size_t size() const { return Nodes.size(); };
 

src/Physics/src/DiscreteRegion.h

@@ -17,6 +17,8 @@ public:
 
     DiscreteRegion(std::shared_ptr<Contour const> region) : Region{region}, Material{Air()} {};
 
+    DiscreteRegion(std::shared_ptr<Contour const> region, std::vector<size_t> elements) : Region{region}, Elements{elements}, Material{Air()} {};
+
     DiscreteRegion(std::shared_ptr<Contour const> region, MaterialProperties material) : Region{region}, Material{material} {};
 
     std::vector<size_t> const &elements() const { return Elements; };

src/Physics/src/FiniteElementMesh.cpp

@@ -2,16 +2,6 @@
 
 template<size_t P, size_t Q>
 FiniteElementMesh<P, Q>::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 ElementOrder > 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)
@@ -28,6 +18,17 @@ FiniteElementMesh<P, Q>::FiniteElementMesh(Mesh &m) {
             Then selection can be implemented by passing the shared_ptr that the user has (similar to how the Mesh operations work)
      */
 
+    // 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 ElementOrder > 1
+    for (size_t i = 0; i != m.size_points(); ++i) {
+        Nodes.emplace_back(m.point(i));
+    }
+
     Boundaries.reserve(m.size_boundary_constraints());
     for (size_t i = 0; i != m.size_boundary_constraints(); ++i) {
         Boundaries.push_back(std::make_shared<DiscreteBoundary>(m.curve(i), m.boundary_nodes(i)));
@@ -44,9 +45,54 @@ FiniteElementMesh<P, Q>::FiniteElementMesh(Mesh &m) {
         Regions[dt.contour()]->push_back(i);
         Triangles.emplace_back(i, m.nodes_from_triangle(dt));
     }
+    */
+
+    // 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
+
+    MeshData md = m.mesh_data(P);
+
+    Nodes.reserve(md.Nodes.size());
+    for (size_t i = 0; i != md.Nodes.size(); ++i) {
+        Nodes.emplace_back(md.Nodes[i]);
+    }
+
+    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));
+    }
+
+    //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());
+    size_t n = 0;
+    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].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();
+    //sort_regions();
 };
 
 template<size_t P, size_t Q>

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::vector<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<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(Mesh &m);
 

src/Physics/src/FiniteElementMeshInterface.cpp

 #include "FiniteElementMeshInterface.h"
 
-FiniteElementMeshInterface::FiniteElementMeshInterface(std::vector<XY> nodes, std::vector<std::shared_ptr<DiscreteRegion>> r, std::vector<std::shared_ptr<DiscreteBoundary>> b)
+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)
         : Nodes(nodes), Regions(r), Boundaries(b) {
     sort_boundaries();
-    sort_regions();
-};
+    //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)
@@ -20,9 +20,10 @@ std::vector<std::shared_ptr<DiscreteBoundary>> FiniteElementMeshInterface::bound
     } else {
         return std::vector<std::shared_ptr<DiscreteBoundary>>();
     }
-};
+}
 
-std::shared_ptr<DiscreteRegion> FiniteElementMeshInterface::region(std::shared_ptr<Contour const> const &c) const {
+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);
@@ -32,6 +33,8 @@ std::shared_ptr<DiscreteRegion> FiniteElementMeshInterface::region(std::shared_p
     } else {
         return std::make_shared<DiscreteRegion>();
     }
+     */
+    return Regions[c];
 }
 
 void FiniteElementMeshInterface::sort_boundaries() {
@@ -39,7 +42,9 @@ void FiniteElementMeshInterface::sort_boundaries() {
     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
+}
+ */
\ No newline at end of file

src/Physics/src/FiniteElementMeshInterface.h

@@ -16,7 +16,7 @@ class FiniteElementMeshInterface {
 public:
     FiniteElementMeshInterface() {};
 
-    FiniteElementMeshInterface(std::vector<XY> nodes, std::vector<std::shared_ptr<DiscreteRegion>> r, std::vector<std::shared_ptr<DiscreteBoundary>> b);
+    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);
 
     virtual ~FiniteElementMeshInterface() {};
 
@@ -38,13 +38,9 @@ public:
 
     auto &regions() { return Regions; };
 
-    auto &region(size_t i) { return Regions[i]; };
-
     auto const &regions() const { return Regions; };
 
-    auto const &region(size_t i) const { return Regions[i]; };
-
-    std::shared_ptr<DiscreteRegion> region(std::shared_ptr<Contour const> const &c) const;
+    std::shared_ptr<DiscreteRegion> region(std::shared_ptr<Contour const> const &c);
 
     DiagonalMatrixGroup diagonal_matrix_group() const { return DiagonalMatrixGroup{quadrature_size(), elements_size()}; };
 
@@ -81,13 +77,17 @@ public:
 protected:
     std::vector<XY> Nodes;
 
-    std::vector<std::shared_ptr<DiscreteRegion>> Regions; // Contains vector of size_t referencing Triangles (and later Quadrilaterals)
+    // TODO: USE SET OR MAP TO STORE REGIONS AND BOUNDARIES BASED ON CONTOUR AND CURVE SHARED_PTRS
+
+    //std::vector<std::shared_ptr<DiscreteRegion>> Regions; // Contains vector of size_t referencing Triangles (and later Quadrilaterals)
+
+    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();
+    // void sort_regions();
 };
 
 #endif //OERSTED_FINITEELEMENTMESHINTERFACE_H

src/Physics/src/Physics.h

@@ -164,7 +164,7 @@ public:
 
             // Calculate field intensity
             for (auto &dr : this->Domain->regions()) {
-                dr->material().h_from_b(dr->elements(), s->Fx, s->Fy, s->dFxdGx, s->dFydGy, s->dFxdGy);
+                dr.second->material().h_from_b(dr.second->elements(), s->Fx, s->Fy, s->dFxdGx, s->dFydGy, s->dFxdGy);
             }
 
             // Calculate residual

src/Physics/src/Triangle.h

@@ -21,6 +21,8 @@ public:
 
     size_t const id() const { return ID; };
 
+    void id(size_t i) { ID = i; };
+
 protected:
     size_t ID;
 };

test/LibraryIntegration/Mesh_To_FEM_Test.cpp

@@ -69,7 +69,8 @@ TEST(Full_Circle, Magnetostatic_Uniform_Current_Density) {
 
     FunctionArguments fargs;
 
-    msph.add_current_density([](FunctionArguments args) { return 1.0 / (2.0 * M_PI * 1e-7); }, {fem->region(0)});
+    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)});
 

test/Physics/Finite_Element_Mesh_Test.cpp

@@ -10,7 +10,7 @@ public:
 
         tri.push_back(Triangle<1>(0, {1, 2, 0}));
 
-        reg.push_back(std::make_shared<DiscreteRegion>(std::vector<size_t>{0}));
+        reg.emplace(c[0], 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}));
@@ -20,7 +20,9 @@ public:
     std::vector<XY> node;
     std::vector<Triangle<1>> tri;
 
-    std::vector<std::shared_ptr<DiscreteRegion>> reg;
+    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;
 
     template<size_t QuadratureOrder>
@@ -81,7 +83,7 @@ TEST_F(MasterTriangle, magnetostatic) {
 
     auto ff = [](FunctionArguments args) { return 1.0 * args["t"]; };
 
-    ms.add_current_density(ff, reg);
+    ms.add_current_density(ff, {reg[c[0]]});
 
     ms.build_matrices();
 
@@ -115,13 +117,15 @@ 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.push_back(std::make_shared<DiscreteRegion>(std::vector<size_t>{0, 1, 2, 3})); // Triangles 4,5,6 are duplicates
+        reg.emplace(c[0], 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<DiscreteRegion>> reg;
+    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;
 
     template<size_t QuadratureOrder>
@@ -244,8 +248,10 @@ public:
         tri.push_back(Triangle<1>(0, {0, 1, 2}));
         tri.push_back(Triangle<1>(1, {1, 3, 2}));
 
-        reg.push_back(std::make_shared<DiscreteRegion>(std::vector<size_t>{0}));
-        reg.push_back(std::make_shared<DiscreteRegion>(std::vector<size_t>{1}));
+        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);
 
         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}));
@@ -262,7 +268,9 @@ public:
     std::vector<XY> node;
     std::vector<Triangle<1>> tri;
 
-    std::vector<std::shared_ptr<DiscreteRegion>> reg;
+    std::vector<std::shared_ptr<const Contour>> c{std::make_shared<const Contour>(), 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;
 
     Eigen::Vector4d vx_m;
@@ -372,7 +380,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[0]});
+    ms.add_current_density(ff, {reg[c[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);
@@ -380,7 +388,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[1]});
+    ms.add_current_density(ff, {reg[c[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);
@@ -388,7 +396,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[0], reg[1]});
+    ms.add_current_density(ff, {reg[c[0]], reg[c[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);
@@ -480,7 +488,7 @@ 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.push_back(inner_region);
+        reg.emplace(c[0], inner_region);
         bnd.push_back(inner_boundary);
 
         tri.push_back(Triangle<1>(6, {1, 7, 2}));
@@ -504,15 +512,20 @@ 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.push_back(outer_region);
+        reg.emplace(c[1], outer_region);
         bnd.push_back(outer_boundary);
 
+        EXPECT_EQ(reg.size(), 2);
+
         femesh = std::make_shared<FiniteElementMesh<1, 1>>(node, tri, reg, bnd);
     }
 
     std::vector<XY> node;
     std::vector<Triangle<1>> tri;
-    std::vector<std::shared_ptr<DiscreteRegion>> reg;
+
+    std::vector<std::shared_ptr<const Contour>> c{std::make_shared<const Contour>(),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::shared_ptr<DiscreteRegion> inner_region;
@@ -658,7 +671,8 @@ TEST_F(TwoRegionHexagon, magnetostatic_forcing_air) {
 }
 
 TEST_F(TwoRegionHexagon, magnetostatic_forcing_1000) {
-    femesh->region(1)->material(Linear1000);
+    //femesh->region(1)->material(Linear1000);
+    outer_region->material(Linear1000);
 
     Magnetostatic ms{femesh};
     ms.build_matrices();
@@ -739,7 +753,7 @@ public:
         tri.push_back(Triangle<1>(1, {0, 2, 3}));
 
         inner_region = std::make_shared<DiscreteRegion>(std::vector<size_t>{0, 1});
-        reg.push_back(inner_region);
+        reg.emplace(c[0], inner_region);
 
         // Outer Region
         r = 2.0;
@@ -757,7 +771,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.push_back(outer_region);
+        reg.emplace(c[1], outer_region);
 
         // Periodic Boundary
         pbi0 = std::make_shared<DiscreteBoundary>(std::vector<size_t>{0,1});
@@ -771,11 +785,16 @@ public:
         bnd.push_back(pbo1);
 
         femesh = std::make_shared<FiniteElementMesh<1, 1>>(node, tri, reg, bnd);
+
+        EXPECT_EQ(reg.size(), 2);
     }
 
     std::vector<XY> node;
     std::vector<Triangle<1>> tri;
-    std::vector<std::shared_ptr<DiscreteRegion>> reg;
+
+    std::vector<std::shared_ptr<const Contour>> c{std::make_shared<const Contour>(),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::shared_ptr<DiscreteRegion> inner_region;
@@ -958,10 +977,10 @@ public:
         tri.push_back(Triangle<1>(5, {0, 6, 1}));
 
         region_ip = std::make_shared<DiscreteRegion>(std::vector<size_t>{0, 1, 2});
-        reg.push_back(region_ip);
+        reg.emplace(c[0], region_ip);
 
         region_in = std::make_shared<DiscreteRegion>(std::vector<size_t>{3, 4, 5});
-        reg.push_back(region_in);
+        reg.emplace(c[1], region_in);
 
         inner_boundary_0 = std::make_shared<DiscreteBoundary