Commit 08f11fbf authored by JasonPries's avatar JasonPries

Add FiniteElementMesh Boundary<2> selection by std::shared_ptr<Curve>

Begin implementation of PeriodicBoundaryCondition construction from PeriodicBoundaryPairs
parent 8d837fbb
......@@ -99,6 +99,7 @@ std::vector<size_t> BoundaryConstraint::nodes(Mesh const &m) {
for (size_t i : DartConstraints) {
n.push_back(m.dart_constraint(i).base(m));
}
n.push_back(m.dart_constraint(DartConstraints.back()).tip(m));
return n;
}
\ No newline at end of file
......@@ -11,8 +11,10 @@ class Boundary {
};
template<>
class Boundary<2> {
class Boundary<2> { // TODO: Rename DiscreteBoundary
public:
Boundary() {};
Boundary(std::vector<size_t> &&nodes) : CurvePtr{nullptr}, Nodes{nodes} {};
Boundary(std::shared_ptr<Curve const> cptr, std::vector<size_t> &&nodes) : CurvePtr{cptr}, Nodes{nodes} {};
......@@ -27,6 +29,8 @@ public:
auto const &node(size_t i) const { return Nodes[i]; };
auto curve() const { return CurvePtr; };
protected:
std::shared_ptr<Curve const> CurvePtr;
std::vector<size_t> Nodes;
......
#include "BoundaryCondition.h"
#include "BoundaryCondition.h"
\ No newline at end of file
......@@ -56,6 +56,32 @@ class PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder> : public Bound
public:
PeriodicBoundaryCondition(std::vector<BoundaryMap> map, bool antiperiodic) : Map{map}, Antiperiodic{antiperiodic} {};
PeriodicBoundaryCondition(std::vector<PeriodicBoundaryPair> const &periodic_boundary_pairs, FiniteElementMesh<2,ElementOrder> const &fem, bool antiperiodic) : Antiperiodic{antiperiodic} {
Map.reserve(periodic_boundary_pairs.size());
double_t sign = (antiperiodic ? -1.0 : 1.0);
for(auto const &pbp : periodic_boundary_pairs) {
// TODO: Change Map type to std::set<VariableMap>, add comparison operator to VariableMap, and make sure node indicies are lexigraphically ordered in VariableMap
auto const &n0 = fem.boundary(pbp.curve0())->nodes();
auto const &n1 = fem.boundary(pbp.curve1())->nodes();
std::vector<VariableMap> vmap;
vmap.reserve(n0.size());
if (pbp.orientation()) {
for (size_t i = 0; i != n0.size(); ++i) {
vmap.emplace_back(n0[i], n1[i], sign);
}
} else {
size_t j = n1.size();
for (size_t i = 0; i != n0.size(); ++i) {
vmap.emplace_back(n0[i], n1[--j], sign);
}
}
Map.emplace_back(vmap);// TODO: needs testing, allows duplicate nodes which may not be handle correctly by constraint algorithm
}
}
void apply(std::vector<Eigen::Triplet<double>> &triplets) const override {
double sign{1.0 - 2.0 * Antiperiodic};
......
......@@ -13,14 +13,18 @@
#include "Region.h"
template<size_t Dimension, size_t Order>
class FiniteElementMesh {};
class FiniteElementMesh {
};
template<size_t Order>
class FiniteElementMesh<2,Order> {
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) {
auto comp = [](std::shared_ptr<Boundary<2>> const &x, std::shared_ptr<Boundary<2>> const &y) { return (size_t) (x->curve().get()) < (size_t) (y->curve().get()); };
std::sort(Boundaries.begin(), Boundaries.end(), comp);
};
FiniteElementMesh(Mesh &m);
......@@ -36,6 +40,19 @@ public:
auto &boundary(size_t i) { return Boundaries[i]; };
std::shared_ptr<Boundary<2>> boundary(std::shared_ptr<Curve const> const &c) const {
auto comp = [](std::shared_ptr<Boundary<2>> const &x, size_t const &y) { return (size_t) (x->curve().get()) < y; };
auto result = std::lower_bound(Boundaries.begin(), Boundaries.end(), (size_t) c.get(), comp);
if (result != Boundaries.end()) {
return *result;
} else {
std::shared_ptr<Boundary<2>> null;
return null;
}
};
auto const &boundaries() const { return Boundaries; };
auto const &nodes() const { return Nodes; };
......@@ -76,7 +93,7 @@ public:
DerivativeMatrixGroup<TriangleQuadrature<Q>::size> derivative() const {
DerivativeMatrixGroup<TriangleQuadrature<Q>::size> df(Nodes.size(), Triangles.size(), Triangle<Order>::NumNodes);
for (size_t i = 0; i!= Triangles.size();++i) {
for (size_t i = 0; i != Triangles.size(); ++i) {
Triangles[i].derivative<Q>(df, Nodes);
}
......@@ -88,7 +105,7 @@ public:
std::vector<std::vector<XY>> qp;
qp.reserve(Triangles.size());
for(size_t i = 0; i!= Triangles.size();++i){
for (size_t i = 0; i != Triangles.size(); ++i) {
qp.push_back(Triangles[i].quadrature_points<Q>(Nodes));
}
......
......@@ -173,6 +173,10 @@ public:
BoundaryConditions.push_back(std::make_shared<PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder>>(map, antiperiodic));
}
void add_periodic_boundary(std::vector<PeriodicBoundaryPair> periodic_boundary_pairs, bool antiperiodic) {
BoundaryConditions.push_back(std::make_shared<PeriodicBoundaryCondition<2, ElementOrder, QuadratureOrder>>(periodic_boundary_pairs, Domain, antiperiodic));
}
template <typename... Args>
auto add_sliding_interface(Args&&... args) {
auto condition = std::make_shared<SlidingInterface<2, ElementOrder, QuadratureOrder>>(std::forward<Args>(args)...);
......
......@@ -7,10 +7,31 @@ class ContinuousBoundaryPair {
public:
ContinuousBoundaryPair(std::shared_ptr<Curve const> c0, std::shared_ptr<Curve const> c1, bool dir) : Curve0{c0}, Curve1{c1}, Orientation{dir} {};
std::shared_ptr<Curve const> curve0() const { return Curve0; };
std::shared_ptr<Curve const> curve1() const { return Curve1; };
bool orientation() const { return Orientation; };
protected:
std::shared_ptr<Curve const> Curve0;
std::shared_ptr<Curve const> Curve1;
bool Orientation;
};
class PeriodicBoundaryPair : public ContinuousBoundaryPair { // TODO: Need to check invariant on construction
public:
PeriodicBoundaryPair(std::shared_ptr<Curve const> c0, std::shared_ptr<Curve const> c1, bool dir, std::shared_ptr<Vertex const> origin, double_t angle) : ContinuousBoundaryPair{c0, c1, dir}, Origin{origin}, Angle{angle} {};
std::shared_ptr<Vertex const> origin() const { return Origin; };
double_t angle() const { return Angle; };
protected:
std::shared_ptr<Vertex const> Origin;
double_t Angle;
};
#endif //OERSTED_CONTINUOUSBOUNDARYPAIR_H
......@@ -203,8 +203,8 @@ void Sketch::save_as<SaveMethod::Rasterize>(std::string path, std::string file_n
fs.close();
}
std::vector<ContinuousBoundaryPair> Sketch::select_periodic_boundary_pairs(std::shared_ptr<Vertex const> const &v0, double_t angle) const {
std::vector<ContinuousBoundaryPair> cbp;
std::vector<PeriodicBoundaryPair> Sketch::select_periodic_boundary_pairs(std::shared_ptr<Vertex const> const &v0, double_t angle) const {
std::vector<PeriodicBoundaryPair> cbp;
std::vector<std::shared_ptr<Curve const>> bc = Boundary->curves();
......@@ -226,10 +226,10 @@ std::vector<ContinuousBoundaryPair> Sketch::select_periodic_boundary_pairs(std::
}
if (m == MatchOrientation::Forward) {
cbp.emplace_back(*c0, *c1, true);
cbp.emplace_back(*c0, *c1, true, v0, angle);
break;
} else if (m == MatchOrientation::Reverse) {
cbp.emplace_back(*c0, *c1, false);
cbp.emplace_back(*c0, *c1, false, v0, angle);
break;
} else {
++c0;
......
......@@ -81,9 +81,9 @@ public:
std::shared_ptr<Vertex const> vertex(size_t i) const { return Verticies[i]; };
std::vector<std::shared_ptr<Curve const>> const curves() { return std::vector<std::shared_ptr<Curve const>>{Curves.begin(),Curves.end()};};
std::vector<std::shared_ptr<Curve const>> const curves() { return std::vector<std::shared_ptr<Curve const>>{Curves.begin(), Curves.end()}; };
std::vector<ContinuousBoundaryPair> select_periodic_boundary_pairs(std::shared_ptr<Vertex const> const &v0, double_t angle) const;
std::vector<PeriodicBoundaryPair> select_periodic_boundary_pairs(std::shared_ptr<Vertex const> const &v0, double_t angle) const;
std::vector<std::shared_ptr<Curve const>> select_radial_boundary(std::shared_ptr<Vertex const> const &v0, double_t radius) const;
......
......@@ -55,9 +55,10 @@ TEST(Library_Integration, Full_Circle_Uniform_Current_Density) {
// Test ordering of boundary nodes
double_t a1{std::atan2(p.y(), p.x())};
if (a1 < 0) {
if (a1 < 0 || (a1 == 0 && a0 > M_PI)) {
a1 += 2 * M_PI;
}
EXPECT_TRUE(a1 > a0);
a0 = a1;
}
......@@ -179,10 +180,10 @@ TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) {
auto periodic_boundary = sk.select_periodic_boundary_pairs(v0, 90.0);
{
for (auto &bp : periodic_boundary) {
if (bp.Curve0.get() == l01.get()) {
EXPECT_EQ(bp.Curve1->is_identical(l01, v0, 90.0), MatchOrientation::Reverse);
} else if (bp.Curve0.get() == l12.get()) {
EXPECT_EQ(bp.Curve1->is_identical(l12, v0, 90.0), MatchOrientation::Reverse);
if (bp.curve0().get() == l01.get()) {
EXPECT_EQ(bp.curve1()->is_identical(l01, v0, 90.0), MatchOrientation::Reverse);
} else if (bp.curve0().get() == l12.get()) {
EXPECT_EQ(bp.curve1()->is_identical(l12, v0, 90.0), MatchOrientation::Reverse);
} else {
GTEST_NONFATAL_FAILURE_("No matching boundary found");
}
......@@ -233,8 +234,56 @@ TEST(Library_Integration, Quadrter_Circle_Mirror_Copy_Uniform_Current_Density) {
me.MinimumElementSize = 0.01;
me.MinimumElementQuality = 0.5;
for (auto &pb : periodic_boundary) {
auto bc0 = me.boundary_constraint(pb.curve0());
if (bc0) {
bc0->uniform_discretization(true);
} else {
GTEST_NONFATAL_FAILURE_("No matching boundary found.");
}
auto bc1 = me.boundary_constraint(pb.curve1());
if (bc1) {
bc1->uniform_discretization(true);
} else {
GTEST_NONFATAL_FAILURE_("No matching boundary found.");
}
}
result = me.refine();
ASSERT_TRUE(result);
me.save_as(SAVE_DIR, std::string("quarter_circle_mirror_copy_mesh"));
// Create FiniteElementMesh
FiniteElementMesh<2,1> fem{me};
for (auto & pb : periodic_boundary) {
auto const &nodes0 = fem.boundary(pb.curve0())->nodes();
auto const &nodes1 = fem.boundary(pb.curve1())->nodes();
for (size_t i = 0; i != nodes0.size(); ++i) {
auto const &n0 = fem.node(nodes0[i]);
auto const &n1 = fem.node(nodes1[nodes0.size() - 1 - i]);
EXPECT_NEAR(n0.y(), 0.0, FLT_EPSILON);
EXPECT_NEAR(n1.x(), 0.0, FLT_EPSILON);
EXPECT_NEAR(n0.x(), n1.y(), FLT_EPSILON);
}
}
// Create 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.add_periodic_boundary(periodic_boundary, false); GTEST_FATAL_FAILURE_("Needs testing, allows duplicate nodes which may not be handle correctly by constraint algorithm");
//msph.build_matrices();
//msph.apply_conditions();
std::cout << std::endl;
}
\ No newline at end of file
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