Commit 1cd258a3 authored by JasonPries's avatar JasonPries
Browse files

Implements FiniteElementMesh<2>::make_discontinuous

Changes a DiscreteBoundary<2> to a DiscontinuousBoundary<2> This will be used for implementing relative motion in a user friendly way
parent 14abfb38
...@@ -154,7 +154,7 @@ public: ...@@ -154,7 +154,7 @@ public:
} }
} }
protected: protected: // TODO: Use VariableMap instead?
std::vector<size_t> First; std::vector<size_t> First;
std::vector<size_t> Second; std::vector<size_t> Second;
std::vector<double> Value; std::vector<double> Value;
......
...@@ -19,9 +19,11 @@ class DiscreteBoundary<2> { ...@@ -19,9 +19,11 @@ class DiscreteBoundary<2> {
public: public:
DiscreteBoundary() {}; DiscreteBoundary() {};
DiscreteBoundary(std::vector<size_t> nodes) : CurvePtr{nullptr}, Nodes{nodes} {}; DiscreteBoundary(std::vector<size_t> nodes, bool closed = false) : CurvePtr{nullptr}, Nodes{nodes}, Closed{closed} {};
DiscreteBoundary(std::shared_ptr<Curve const> cptr, std::vector<size_t> nodes) : CurvePtr{cptr}, Nodes{nodes} {}; DiscreteBoundary(std::shared_ptr<Curve const> cptr, std::vector<size_t> nodes, bool closed = false) : CurvePtr{cptr}, Nodes{nodes}, Closed{closed} {};
virtual ~DiscreteBoundary() {};
size_t size() const { return Nodes.size(); }; size_t size() const { return Nodes.size(); };
...@@ -35,9 +37,12 @@ public: ...@@ -35,9 +37,12 @@ public:
std::shared_ptr<Curve const> curve() const { return CurvePtr; }; std::shared_ptr<Curve const> curve() const { return CurvePtr; };
bool is_closed() const { return Closed; };
protected: protected:
std::shared_ptr<Curve const> CurvePtr; std::shared_ptr<Curve const> CurvePtr;
std::vector<size_t> Nodes; std::vector<size_t> Nodes;
bool Closed;
}; };
template<> template<>
...@@ -45,9 +50,9 @@ class DiscontinuousBoundary<2> : public DiscreteBoundary<2> { ...@@ -45,9 +50,9 @@ class DiscontinuousBoundary<2> : public DiscreteBoundary<2> {
public: public:
DiscontinuousBoundary() {}; DiscontinuousBoundary() {};
DiscontinuousBoundary(std::vector<size_t> nodes, std::vector<size_t> dnodes, std::vector<double_t> vals) : DiscreteBoundary{nodes}, DiscontinuousNodes{dnodes}, Values{vals} {}; DiscontinuousBoundary(std::vector<size_t> nodes, std::vector<size_t> dnodes, std::vector<double_t> vals, bool closed = false) : DiscreteBoundary{nodes, closed}, DiscontinuousNodes{dnodes}, Values{vals} {};
DiscontinuousBoundary(std::shared_ptr<Curve const> cptr, std::vector<size_t> nodes, std::vector<size_t> dnodes, std::vector<double_t> vals) : DiscreteBoundary{cptr, nodes}, DiscontinuousNodes{dnodes}, Values{vals} {}; DiscontinuousBoundary(std::shared_ptr<Curve const> cptr, std::vector<size_t> nodes, std::vector<size_t> dnodes, std::vector<double_t> vals, bool closed = false) : DiscreteBoundary{cptr, nodes, closed}, DiscontinuousNodes{dnodes}, Values{vals} {};
std::vector<size_t> dnodes() { return DiscontinuousNodes; }; std::vector<size_t> dnodes() { return DiscontinuousNodes; };
...@@ -59,7 +64,7 @@ public: ...@@ -59,7 +64,7 @@ public:
size_t dnode(size_t i) { return DiscontinuousNodes[i]; }; size_t dnode(size_t i) { return DiscontinuousNodes[i]; };
size_t const &dnode(size_t i) const {return DiscontinuousNodes[i]; }; size_t const &dnode(size_t i) const { return DiscontinuousNodes[i]; };
protected: protected:
std::vector<size_t> DiscontinuousNodes; std::vector<size_t> DiscontinuousNodes;
...@@ -87,7 +92,7 @@ protected: ...@@ -87,7 +92,7 @@ protected:
class VariableMap { class VariableMap {
public: public:
VariableMap(size_t first, size_t second, double value) : First{std::min(first,second)}, Second{std::max(first,second)}, Value{value} {}; VariableMap(size_t first, size_t second, double value) : First{std::min(first, second)}, Second{std::max(first, second)}, Value{value} {};
size_t first() const { return First; }; size_t first() const { return First; };
......
...@@ -53,6 +53,77 @@ public: ...@@ -53,6 +53,77 @@ public:
} }
}; };
std::shared_ptr<DiscontinuousBoundary<2>> make_discontinuous(std::shared_ptr<Curve const> const &c) {
// TODO: Should this be implemented in the refineable mesh instead? It might be simpler.
std::shared_ptr<DiscreteBoundary<2>> b = boundary(c);
if (!b) {
std::shared_ptr<DiscontinuousBoundary<2>> null;
return null;
}
std::shared_ptr<DiscontinuousBoundary<2>> db = std::dynamic_pointer_cast<DiscontinuousBoundary<2>>(b);
if (db) {
return db;
}
// Copy nodes
std::vector<size_t> bnodes = b->nodes();
std::vector<size_t> dnodes;
dnodes.reserve(bnodes.size());
for (size_t i = 0; i != bnodes.size(); ++i) {
dnodes.push_back(Nodes.size());
Nodes.emplace_back(Nodes[bnodes[i]]);
}
// Construct new DiscontinuousBoundary and replace
db = std::make_shared<DiscontinuousBoundary<2>>(bnodes, dnodes, std::vector<double_t>(dnodes.size(), 1.0), b->is_closed());
std::replace(Boundaries.begin(), Boundaries.end(), b, std::dynamic_pointer_cast<DiscreteBoundary<2>>(db));
// Renumber nodes in elements when the contains nodes with the same orientation as the boundary
// TODO: This a brute renumbering algorithm with O(t*n) complexity where t = Triangles.size() and n = bnodes.size()
// TODO: Sort + renumber should result in O(n*log(t)) complexity, but all present assumptions about this class assume the the ordering of the elements never changes
for (size_t i = 0; i != (dnodes.size() - 1 + db->is_closed()); ++i) {
size_t ii = (i + 1) % dnodes.size();
double_t dx0 = Nodes[bnodes[ii]].x() - Nodes[bnodes[i]].x();
double_t dy0 = Nodes[bnodes[ii]].y() - Nodes[bnodes[i]].y();
for (Triangle<Order> &t : Triangles) {
for (size_t k = 0; k != 3; ++k) {
if (t.node(k) == bnodes[i]) {
size_t k0 = t.node(0);
size_t k1 = t.node(1);
size_t k2 = t.node(2);
double_t dx1 = (Nodes[k0].x() + Nodes[k1].x() + Nodes[k2].x()) / 3.0 - (Nodes[bnodes[i]].x() + Nodes[bnodes[ii]].x()) / 2.0;
double_t dy1 = (Nodes[k0].y() + Nodes[k1].y() + Nodes[k2].y()) / 3.0 - (Nodes[bnodes[i]].y() + Nodes[bnodes[ii]].y()) / 2.0;
double_t cross = dx0 * dy1 - dy0 * dx1;
if (cross < 0) {
t.node(k, dnodes[i]);
}
}
if (t.node(k) == bnodes[ii]) {
size_t k0 = t.node(0);
size_t k1 = t.node(1);
size_t k2 = t.node(2);
double_t dx1 = (Nodes[k0].x() + Nodes[k1].x() + Nodes[k2].x()) / 3.0 - (Nodes[bnodes[i]].x() + Nodes[bnodes[ii]].x()) / 2.0;
double_t dy1 = (Nodes[k0].y() + Nodes[k1].y() + Nodes[k2].y()) / 3.0 - (Nodes[bnodes[i]].y() + Nodes[bnodes[ii]].y()) / 2.0;
double_t cross = dx0 * dy1 - dy0 * dx1;
if (cross < 0) {
t.node(k, dnodes[ii]);
}
}
}
}
}
return db;
}
auto const &boundaries() const { return Boundaries; }; auto const &boundaries() const { return Boundaries; };
auto const &nodes() const { return Nodes; }; auto const &nodes() const { return Nodes; };
...@@ -112,7 +183,6 @@ public: ...@@ -112,7 +183,6 @@ public:
return qp; return qp;
} }
//void save_as(std::string path, std::string file_name) const;
void write_scalar(Eigen::VectorXd const &v, std::string path, std::string file_name) const { void write_scalar(Eigen::VectorXd const &v, std::string path, std::string file_name) const {
if (!boost::filesystem::exists(path)) { if (!boost::filesystem::exists(path)) {
boost::filesystem::create_directories(path); boost::filesystem::create_directories(path);
...@@ -145,19 +215,10 @@ public: ...@@ -145,19 +215,10 @@ public:
} }
// Write values // Write values
for (size_t i = 0; i!=v.size(); ++i){ for (size_t i = 0; i != v.size(); ++i) {
fs << v(i) << std::endl; fs << v(i) << std::endl;
} }
/*
for (size_t e = 0; e != Edges.size(); ++e) {
Point const v0 = base(e);
Point const v1 = base(next(e));
Point const v2 = base(next(next(e)));
fs << v0.X << ',' << v1.X << ',' << v2.X << ',' << v0.Y << ',' << v1.Y << ',' << v2.Y << '\n';
}
*/
fs.close(); fs.close();
} }
......
...@@ -26,10 +26,12 @@ class SmallVector<2> { ...@@ -26,10 +26,12 @@ class SmallVector<2> {
public: public:
SmallVector() : Value{} {}; SmallVector() : Value{} {};
SmallVector(Point const& p) : Value{p.X, p.Y} {};
SmallVector(double x, double y) : Value{x, y} {}; SmallVector(double x, double y) : Value{x, y} {};
SmallVector(SmallVector<2> const &v) : Value{v.x(), v.y()} {};
SmallVector(Point const& p) : Value{p.X, p.Y} {};
double value(size_t i) const { return Value[i]; }; double value(size_t i) const { return Value[i]; };
double x() const { return Value[0]; }; double x() const { return Value[0]; };
...@@ -53,6 +55,8 @@ public: ...@@ -53,6 +55,8 @@ public:
SmallVector(double x, double y, double z) : Value{x, y, z} {}; SmallVector(double x, double y, double z) : Value{x, y, z} {};
SmallVector(SmallVector<3> const &v) : Value{v.x(), v.y(), v.z()} {};
double value(size_t i) const { return Value[i]; } double value(size_t i) const { return Value[i]; }
double x() const { return Value[0]; } double x() const { return Value[0]; }
......
...@@ -194,11 +194,6 @@ public: ...@@ -194,11 +194,6 @@ public:
* Instead, accept a std::shared_ptr<Curve> * Instead, accept a std::shared_ptr<Curve>
* Call Domain.make_discontinuous(std::shared_ptr<Curve>) * Call Domain.make_discontinuous(std::shared_ptr<Curve>)
* Construct sliding_interface with the returned std::shared_ptr<DiscontinuousBoundary> * Construct sliding_interface with the returned std::shared_ptr<DiscontinuousBoundary>
*
* TODO: Implement Domain.make_discontinuous(std::shared_ptr<Curve>)
* Try casting the std::shared_ptr<DiscreteBoundary> to std::shared_ptr<DiscontinuousBoundary>
* If successful, return the std::shared_ptr<DiscontinuousBoundary>
* Else, construct DiscontinuousBoundary, replace the existing std::shared_ptr<DiscreteBoundary>, and return the new std::shared_ptr<DiscontinuousBoundary>
*/ */
auto condition = std::make_shared<SlidingInterface<2, ElementOrder, QuadratureOrder>>(std::forward<Args>(args)...); auto condition = std::make_shared<SlidingInterface<2, ElementOrder, QuadratureOrder>>(std::forward<Args>(args)...);
BoundaryConditions.push_back(condition); BoundaryConditions.push_back(condition);
......
...@@ -44,6 +44,8 @@ public: ...@@ -44,6 +44,8 @@ public:
} }
}; };
void node(size_t const &i, size_t n) { Node[i] = n; };
size_t const &node(size_t const &i) const { return Node[i]; }; size_t const &node(size_t const &i) const { return Node[i]; };
template<size_t D> template<size_t D>
......
...@@ -1118,49 +1118,47 @@ public: ...@@ -1118,49 +1118,47 @@ public:
region_ip = std::make_shared<Region<2>>(std::vector<size_t>{0, 1, 2}); 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}); region_in = std::make_shared<Region<2>>(std::vector<size_t>{3, 4, 5});
std::shared_ptr<CircularArc const> ci{};
boundary_i = std::make_shared<DiscreteBoundary<2>>(ci, std::vector<size_t>{1, 2, 3, 4, 5, 6}, true);
reg.push_back(region_ip); reg.push_back(region_ip);
reg.push_back(region_in); reg.push_back(region_in);
bnd.push_back(boundary_i);
// Outer Region // Outer Region
r = 1.0;
for (size_t i = 0; i != 6; ++i) {
double a = (i * M_PI) / 3.0;
node.push_back(XY{r * std::cos(a), r * std::sin(a)});
}
r = 2.0; r = 2.0;
for (size_t i = 0; i != 6; ++i) { for (size_t i = 0; i != 6; ++i) {
double a = (i * M_PI) / 3.0 + M_PI / 6.0; double a = (i * M_PI) / 3.0 + M_PI / 6.0;
node.push_back(XY{r * std::cos(a), r * std::sin(a)}); node.push_back(XY{r * std::cos(a), r * std::sin(a)});
} }
tri.push_back(Triangle<1>(6, {7, 13, 8})); tri.push_back(Triangle<1>(6, {1, 7, 2}));
tri.push_back(Triangle<1>(7, {8, 13, 14})); tri.push_back(Triangle<1>(7, {2, 7, 8}));
tri.push_back(Triangle<1>(8, {8, 14, 9})); tri.push_back(Triangle<1>(8, {2, 8, 3}));
tri.push_back(Triangle<1>(9, {9, 14, 15})); tri.push_back(Triangle<1>(9, {3, 8, 9}));
tri.push_back(Triangle<1>(10, {9, 15, 10})); tri.push_back(Triangle<1>(10, {3, 9, 4}));
tri.push_back(Triangle<1>(11, {10, 15, 16})); tri.push_back(Triangle<1>(11, {4, 9, 10}));
tri.push_back(Triangle<1>(12, {10, 16, 11})); tri.push_back(Triangle<1>(12, {4, 10, 5}));
tri.push_back(Triangle<1>(13, {11, 16, 17})); tri.push_back(Triangle<1>(13, {5, 10, 11}));
tri.push_back(Triangle<1>(14, {11, 17, 12})); tri.push_back(Triangle<1>(14, {5, 11, 6}));
tri.push_back(Triangle<1>(15, {12, 17, 18})); tri.push_back(Triangle<1>(15, {6, 11, 12}));
tri.push_back(Triangle<1>(16, {12, 18, 7})); tri.push_back(Triangle<1>(16, {6, 12, 1}));
tri.push_back(Triangle<1>(17, {7, 18, 13})); tri.push_back(Triangle<1>(17, {1, 12, 7}));
region_o = std::make_shared<Region<2>>(std::vector<size_t>{6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17}); region_o = std::make_shared<Region<2>>(std::vector<size_t>{6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17});
boundary_oo = std::make_shared<DiscreteBoundary<2>>(std::vector<size_t>{13, 14, 15, 16, 17, 18}); boundary_o = std::make_shared<DiscreteBoundary<2>>(std::vector<size_t>{7, 8, 9, 10, 11, 12}, false);
sliding_boundary = std::make_shared<DiscontinuousBoundary<2>>(std::vector<size_t>{1, 2, 3, 4, 5, 6}, std::vector<size_t>{7, 8, 9, 10, 11, 12}, std::vector<double>(6, 1.0));
reg.push_back(region_o); reg.push_back(region_o);
bnd.push_back(boundary_oo); bnd.push_back(boundary_o);
bnd.push_back(sliding_boundary);
femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd); femesh = FiniteElementMesh<2, 1>(node, tri, reg, bnd);
sliding_boundary = femesh.make_discontinuous(ci);
} }
std::vector<XY> node; std::vector<XY> node;
...@@ -1173,7 +1171,8 @@ public: ...@@ -1173,7 +1171,8 @@ public:
std::shared_ptr<Region<2>> region_o; std::shared_ptr<Region<2>> region_o;
std::shared_ptr<DiscontinuousBoundary<2>> sliding_boundary; std::shared_ptr<DiscontinuousBoundary<2>> sliding_boundary;
std::shared_ptr<DiscreteBoundary<2>> boundary_oo; std::shared_ptr<DiscreteBoundary<2>> boundary_o;
std::shared_ptr<DiscreteBoundary<2>> boundary_i;
FiniteElementMesh<2, 1> femesh; FiniteElementMesh<2, 1> femesh;
...@@ -1246,7 +1245,7 @@ TEST_F(TwoRegionHexagonRotation, magnetostatic_forcing) { ...@@ -1246,7 +1245,7 @@ TEST_F(TwoRegionHexagonRotation, magnetostatic_forcing) {
ms.add_current_density(f_positive, {region_ip}); ms.add_current_density(f_positive, {region_ip});
ms.add_current_density(f_negative, {region_in}); ms.add_current_density(f_negative, {region_in});
ms.add_magnetic_insulation({boundary_oo}); ms.add_magnetic_insulation({boundary_o});
auto position = ms.add_sliding_interface(sliding_boundary, false); auto position = ms.add_sliding_interface(sliding_boundary, false);
std::vector<double> v_expected{0.0, 0.0, 1.15908e-7, 1.15908e-7, 0.0, -1.15908e-7, -1.15908e-7}; std::vector<double> v_expected{0.0, 0.0, 1.15908e-7, 1.15908e-7, 0.0, -1.15908e-7, -1.15908e-7};
...@@ -1288,6 +1287,9 @@ TEST_F(TwoRegionHexagonRotation, magnetostatic_forcing) { ...@@ -1288,6 +1287,9 @@ TEST_F(TwoRegionHexagonRotation, magnetostatic_forcing) {
EXPECT_NEAR(v(i), v_expected[i], FLT_EPSILON); EXPECT_NEAR(v(i), v_expected[i], FLT_EPSILON);
} }
Eigen::VectorXd vv = ms.recover_boundary(v);
femesh.write_scalar(vv, SAVE_DIR, std::string("TwoRegionHexagonRotation") + std::to_string(iter));
// Test flux-density values // Test flux-density values
Eigen::ArrayXd Bx = ms.derivatives().dy(0).transpose() * v; Eigen::ArrayXd Bx = ms.derivatives().dy(0).transpose() * v;
Eigen::ArrayXd By = -ms.derivatives().dx(0).transpose() * v; Eigen::ArrayXd By = -ms.derivatives().dx(0).transpose() * v;
......
...@@ -7,4 +7,6 @@ ...@@ -7,4 +7,6 @@
#include "Physics.hpp" #include "Physics.hpp"
#define SAVE_DIR "./test/output/Physics/"
#endif //OERSTED_TEST_PHYSICS_HPP #endif //OERSTED_TEST_PHYSICS_HPP
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment